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Feedback spring so that EHSV flow is proportional to FADEC 
commanded torque motor current 
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way” valves and discrete (on/off) valves. 
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Summary of HMUX system testing : 
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Summary of HMUX system testing (cont.) : 
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HMUX System Rig Test at WAES 
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P Aircraft Engines Advanced Control Systems 

HMUX FMV Performance: 
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Performance During Throttle Chop 
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the four fuel zones are progressively uncovered. There is sufficient stroke between stages to ensure that ; 
the valves will be staged correctly when the master valve is positioned at the nominal stroke. Two 
dynamic seals are utilized to prevent fuel leakage into the unstaged zones. If an unstaged zone gets 
leakage, fuel coking and high UHC emissions will result. A series of endurance tests were performed on 
various seal designs to select one with the best combination of low repeatable friction and low leakage. 
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NASA Critical Propulsion and Noise Technologies 


The Controlled Pressure Fuel Nozzle System (CPFN) system designed, fabricated and 
demonstrated under the AST program has proven the viability of utilizing CPFN technology in the control 
systems of complex low emission combustors. The comprehensive component-level testing performed at 
Parker, and system integration testing performed at GE using a system rig running closed loop with an 
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The components developed for this program were prototypes and utilized machinings, rather than castings, 
for the housing. The eventual production design will utilize a cast housing and will be much more compac 
and lighter in weight. As the picture shows, the first stage assembly consisting of the force motor, rotating 
and fixed sleeves, and fuel turbine is located in the center of the assembly. The 15 second stage throttling 
valves are placed around the perimeter of the assembly. 
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Typical Step Response 
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The picture above shows the FMUX system being tested with cold ft 
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Abstract 


In order to enable higher performing feedback control and to better evaluate engine quality, 
it is desirable to be able to estimate both the state vector and key engine parameter variations 
and input biases throughout the flight envelope. This is a difficult problem in practice because 
there are less sensors than parameter variations and input biases, some of the sensors are very 
noisy, and some of the sensor measurements are biased. As shown in this report, due to the first 
difficulty, this estimation problem requires the simultaneous use of data from several operating 
points. Hence, a “stacked” system model, that can be used to incorporate data from several 
operating points, is proposed. The filtering is accomplished by using both the Kalman filter 
and its algebraic equivalent, the Information filter. As shown in simulation results, the Kalman 
filter sometimes diverged due to an ill-conditioned covariance matrix (i.e. the elements of the 
covariance matrix grew without bound) whereas the Information filter always converged. It 
was demonstrated that using the “stacked” configuration with data from three, appropriately 
chosen, operating points (from a list of 59 operating points) 15 parameter variations and input 
biases can be estimated accurately. The use of data from four operating points does not allow 
more parameter variations and input biases to be estimated while these 15 parameters variations 
and input biases cannot be estimated accurately with data from one or two operating points. 
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1. Introduction 


Model-based control is expected to provide a significant improvement in the control per- 
formance for General Electric jet engines. One of the key uses of an aircraft engine model 
is online state estimation as a means of estimating the unmeasurable performance outputs of 
the engine such as thrust and stall margin. These estimates can then be used as the current 
values of these performance parameters in a gain-scheduled, servo control scheme. This type 
of control scheme can substantially decrease the conservatism inherent in current approaches 
that directly control measurable output variables which are only indirectly related to the real 
performance outputs. 

In order to determine engine quality it is also desired to estimate certain engine parameter 
and input biases. (Note that “parameter bias” here refers to a parameter variation.) To accom- 
plish both this objective and the objective discussed above, it is necessary to simultaneously 
estimate the state vector x and the vector of engine parameter and input biases which will be 
denoted by u b . Note that some of the elements of u b may actually vary at different operating 
points, but in this study we will assume that these parameters are constant over all operating 
conditions. Even then, the estimation problem is nontrivial because there are fewer elements 
of the sensor measurement vector y than elements of the vector of parameter and input biases 
ii h. some of the sensors are very noisy, and some of the sensor measurements are biased. 

In this study it is initially assumed that we are interested in using steady-state data at 
various operating points to estimate both the steady-state, state vector x„ and u b . The es- 
timation problem at one (generic) operating point is studied under the ideal conditions of no 
plant disturbance and no plant noise. It will be seen that this analysis indicates that data 
from one operating point is not sufficient to uniquely determine x, s and u h . Hence, the use 
of data from more than one operating point is considered. Guidelines are developed on the 
minimum number of operating points needed and a rank condition on a certain matrix is used 
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to determine appropriate combinations of operating points. A “stacked” model structure that 
utilizes data from several operating points is then proposed as the basis for the filter design. 

A key issue in filtering is numerical robustness. It was observed that when the stacked 
structure was used the covariance matrix would occasionally become divergent (i.e. some of 
the elements grew without bound) causing the Kalman filter to diverge. Hence the Information 
(or inverse covariance) filter was used to increase the numerical robustness of the estimation 
process. To familiarize the reader with the Information filter it is described and compared with 
the Kalman filter. 

This report is organized as follows. Section 2 discusses the steady-state estimation problem 

*• 

under the ideal conditions of no plant disturbances or sensor noise. It then proposes a filtering 
structure that uses data at several operating points simultaneously. Section 3 describes state 
and parameter estimation using the Kalman and Information filters. Section 4 illustrates the 
results with a simplified estimation problem. Section 5 presents results for a realistic set of 
parameters and sensors. Finally, Section 6 presents conclusions, including a brief discussion of 
how the results may be extended to estimation using transient data. 

2. Fundamental Problem and Proposed Solution 

To estimate the states and parameters at a particular operating point a linear model has 

to be constructed. Linear models with any combinations of outputs and inputs desired, are 

constructed by linearizing the nonlinear model provided by GE about specific operating points 
and nominal parameter values. If plant disturbances and sensor noise are ignored, then the 
resulting linear model (after discretization) is of the form 

x(k + 1) = Ax{k) + Bu{k) + B b u b (2.1) 

z(k ) = Cx(k) + Du(k) + D b u b (2.2) 

where u denotes the system inputs and u b denotes the parameter variations and input biases. 
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To consider the steady-state behavior of the system (2.1)-(2.2), define 


x 3s = lim x(k), z„ = lim z{k), u„ = lim u(k). (2.3) 

k -400 k—*oo k—>oc 

Then, {2.1 )-(2.2) yield the steady-state equations 


or equivalently 


where 


x as = Ax„ + Bu s , 4- B b u b 
z sa - Cx, s 4- Du,, 4- D b u b 

AX — Z 


A = 


(A -I) B b 
C D b 


u b 


Bu aa 

z aa Du aa 


(2.4) 

(2.5) 

(2.6) 

(2.7) 


Let us assume that x € 1R 3 , z € IR 8 , u b € IR 22 such that X € IR 25 , Z € IR 11 , and 
A € IR U>:25 . (These dimensions are equal to the data given for the GE jet engine considered 
in t his study.) Notice that if (2.6) has a solution, it is not unique. In fact, if a solution exists, 
then there are an infinite number of solutions. If a solution exist, a necessary condition for a 
unique solution is the satisfaction of the following condition. 

Condition 1. (A necessary condition.) nrcrw(A) > ncol(A) where nrow denotes the number 
of rows and ncol denotes the number of columns. 


A necessary and sufficient condition is the satisfaction of the following. 

Condition 2. ( A necessary and sufficient condition.) rank(A) = ncol(A) where rank denotes 
the rank of the matrix. (This condition requires A to be full column rank.) 

In theory it is possible to satisfy Condition 1 by adding more sensor outputs. For example, 
assume that as before x € IR 3 and u b 6 IR 22 but now z g IR 22 . Then, X € IR 25 , Z € IR 25 , and 
A € IR 25x25 . Since A is square. Condition 1 is now satisfied. Condition 2 is satisfied (such that 
X is unique) only if A is also invertible. Unfortunately, due to engineering limitations in the jet 
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engine problem considered here, it is not possible to sufficiently increase the number of sensors 
such. that dim(y) > dim(ub) where dim denotes vector dimension. Even if this were possible, 
Condition 2 may not be satisfied. Because in practice dim(y) < dim(ub), it is necessary to use 
data from several operating points to perform the estimation. 

Consider N operating points and assume that the parameter and input bias vector Ub does 
not vary with respect to the operating point. Let (Ai, Bi,Q, Di) denote the linear system and 
.r, i/j and Uj, S3 respectively denote the steady-state, state, output, and input vectors at the 
i th operating point. Then, the steady-state conditions at the N operating points are described 
by the following equations. 

« 

•t-l.j* = d" Billies + Bb'iUb'i 

z\ ,« = C\X\' Si + D\U\', a + Db,\Ub,i, 

%2,S3 — A2X2, „ + B2U2,33 + Db % 2 u b,2 
z 2,33 = C 2 X2,,3 + D2U2.33 + A>,2U&,2i 


ZN.ss — AyXy^ aa 4“ By XI y <33 “h Db.S'U-b.S 
z N, 33 ~ CyXy : , 3 + DyUy' S , 4 - Db,yUb,y- 

Concatenating the above equations yields the linear system 


AX = Z 


( 2 . 8 ) 


whore 


Bb, 1 

{A, -I) 

0 

0 

A>,1 

c, 

0 

0 

Bb, 2 

0 

(A 2 - 1) 

0 

Db,2 

0 

C 2 

0 

B b ,y 

0 

0 

••• {At, -I) 

. Db.y 

0 

0 

Cy 


(2.9) 
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-Biu hs , 

u b 


z i,t» ~ D\U\ >a , 

•El ,13 


-BiU2,ts 

X 2 ,as 

, 2 = 

*2,s a ~ A2«2,i» 

* 

•0 

■■ Z 
H 

i 


—BnUn,33 
. V,« ~ DnU2,, s . 


Note that if u b € IR 22 and x< £ IR 3 , Zi £ 1R 8 for i £ {1, 2, JV}, then X £ IR 22+3 " , 2 € IR 11V . 
and A € IR 11Nx(22+3N) . As seen by Table 1, for N > 3 Condition 1 is satisfied. Condition 


2 must be satisfied by choosing an appropriate combination of operating conditions. These 
mrphiTiatinns are dependent upon the physics of the problem and will have to be determined 
experimentally. In practice, given a set of N operating conditions where N is sufficiently lafge 
to satisfy Condition 1, it is necessary to actually construct the matrix A defined in (2.9) to 


determine whether Condition 2 is satisfied. 

If no noise is present, it would be possible to solve for the desired steady-state, state vectors 
x, t , and unknown bias vector u b by simply solving the linear system (2.8). However, since 
noise is present in the real problem it is necessary to use a filtering technique such as Kalman 
or Information filtering. If N operating points are chosen, then the filtering should be based 
on the following “stacked” linear system structure. 


Xi(k + 1) 


‘ A ! 

0 ••• 

0 ■ 


' Xi(fc) ' 


’ By ' 


' B b ,y ' 

x 2 (fc + 1) 

= 

0 

a 2 

0 


x 2 {k) 

+ 

b 2 

u(k) + 

B b , 2 

. £tf(fc + 1) . 


. 0 

... o 

An . 


. x N {k) . 


. Bn . 


. B b ,N . 


' Zi{k) ' 


Cy 

0 ••• 

0 ' 


' Xy(k) ' 


' Dy ■ 


A>, i 



= 

0 

c 2 ••• 

0 


x 2 {k) 

+ 

d 2 

u(k) 4- 

D b , 2 


. Zft(k) 


. 0 

... o 

Cn . 


. x N (k) . 


.Dn . 


. B>b,N . 


u b (2.11) 


u b . 


( 2 . 12 ) 


In practice this filter structure will require the storage of data at several operating points. 
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DE 

mow(A) 

ncol(.4) 

IB 

11 

25 

IB 

22 

28 

111 

33 

31 

111 

44 

34 

IB 

55 

37 

IB 

66 

40 

IB 

77 

43 

IB 

88 

46 

LL 

99 

49 

im 

110 

52 

LiL 

121 

55 

m 

132 

58 


Table 1: Number of Rows and Columns of A 


3. State and Parameter Estimation 

This section begins by showing that the problem of simultaneous input and state estimation 
can be reformulated as a standard state estimation problem. State estimation by standard 
Kalman filtering is then briefly described. Next, state estimation by Information filtering is 
discussed and then compared with Kalman filtering. 

3.1. Input and State Estimation 

It is desired to use Kalman or Information filtering to estimate the state vector x and the 
vector of parameter and input biases u* using a model of the form 

x(k + 1) = Ax(k) + Bu{k ) + J3&t4 (3.1) 

z(fc) = Cx(k) + Du{k) + DtUt. (3.2) 
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This may be approached by expressing (3.1)-(3.2) in a form in which u b is part of an augmented 
state vector [5]. In particular, since u b {k + 1) = u b (k), (3.1)-(3.2) may be expressed as 


"l" 1) Aa U gX aU g(k) + B aU gU(k) 

(3.3) 

z(k) — C aU gX mi g(k) + D au gu(k) 

(3.4) 


where 

Xaug(k) — 


x{k) 

u b 




A B b 
0 I 


Bang — 


B 

0 


, C au g — ^ C D b j . (3.5) 


Hence, in what follows, without loss of generality, only state estimation will be considered. 

3.2. The Kalman Filter 


Consider the time-varying, discrete-time, linear stochastic system 

x(k -f 1) = F(k)x(k) + G(k)u(k ) + w(k) (3.6) 

z(k) = H(k)x(k) + J{k)u(k ) + v(k) (3.7) 

where w(k) is white process noise with covariance Q(k) and v(k) is white sensor noise with 
covariance R{k). 

The Kalman filter is a recursive linear estimator which successively calculates a minimum 
variance estimate for the state x given the observations z. The Kalman filter minimizes the 
mean square estimation error and is optimal with respect to a variety of important criteria 
under specific assumptions about process and observation noise. The Kalman filter can be 
summarized in the following two stages. 

Prediction 


x(k\k - 1) = F(k)x{k - 1 | k - 1) + G{k)u{k) (3.8) 

P(k\k - 1) = F(k)P(k - 1|* - l)F T (k) + Q(k) (3.9) 
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Estimation 


*(*!*) = [/ - (*)]x(*|A: - 1) + W(k)[z(k) - J(k)u{k)} (3.10) 

P(k\k) = P(k\k - 1) - W(k)S(k)W T (k) (3.11) 

where the gain matrix is given by 

W(k) = P(k\k - l)^*)^*) (3.12) 

and the innovation covariance by 

S(k) = H(k)P(k\k - l)H T (k) + R(k). (3.1?) 


Prom (3.10), the Kalman filter state estimate can be interpreted as a linear weighted sum of 
the state prediction and observation. The weights in this averaging process are I - W(k)H(k) 
associated with the prediction and W(fc) associated with the observation. The values of the 
weights depend on the balance of confidence in prediction and observation as specified by the 
process and observation noise covariances. 

The Constant Gain K alman Filter 

For the jet engine problem considered here (and many other problems) the system matrices 
F, G, H, K are constant. In addition, it is sensical to assume the noise covariances Q and R 
constant. In this case, to save processor throughout, it is advantageous to use the steady-state 
value of the Kalman gain W in the filter implementation. The Kalman filter algorithm then 
reduces to the following. 


x(k\k — 1) = Fx(k — l|fc — 1) + Gu(k) 


(3.14) 


x{k\k) = [/ - WH]x(k\k - 1) + W[z(k) - Ju{k )] (3.15) 
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where 


W = PH T S~ l (3.16) 

S = HPH t + R (3.17) 

and P is the nonnegative definite solution of the algebraic Riccati equation, 

P = FPF T + P- FPH t {R + HPH t )~ 1 HPF t . (3.18) 

3.3. The Information Filter 

The Information filter is essentially a Kalman filter expressed in terms of measures of in- 
formation about the parameters (states) of interest rather than direct state estimates and 
their associated covariances. This filter has also been called the inverse covariance form of 
the Kalman filter [3,6, 12]. In this section, the contextual meaning of the information filter is 
derived. 

' For the Information filter there are two key information-analytic variables required, the 
information matrix and information state vector. The information matrix is defined as the 
inverse of the covariance matrix, 

rm = P'HHj).) (3.19) 

The information state vector is a product of the inverse of the covariance matrix (information 
matrix) and the state estimate. 

V(*ljf) = F -1 (*|j)x(i|,7). (3.20) 

The variables, Y{i\j) and y(i\j), form the basis of the information space ideas that are central 
to the estimation presented in this report. For a detail explanation of the derivation of the 
information filter refer to [7]. The Information filter algorithm is summarized as follows: 
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Prediction 


y(k\k - 1) = L(k\k - 1 )y(k - Ij* - 1) (3.21) 

V(k!k - 1) = [F(Ar)y- l (i - 1|* - 1 )F T (k) + Q(Jt)]- 1 (3.22) 

Estimation 

$(*l*) = »(*|* - 1) + H T (k)Rr\k)[z(k) - Ju(Jk)] (3.23) 

Y(k)(k) = Y(k\k - 1) + H T (k)R- l (k)H(k) (3.24) 

where the propagation coefficient is given by 

L(k\k - 1) = F(fc|* - VFWY-'ik - 1|* - 1). (3.2b) 


This is the information form of the Kalman filter, [4,8]. Despite its potential applications, it 
is not widely used and is thinly covered in literature. Barshalom [3] and Maybeck [6] briefly 
discuss the idea of information estimation, but do not explicitly derive the algorithm in terms 
of information as done above, nor do they use it as a principal filtering method. 

Constant Coefficient Information Filter 

As for constant gain Kalman filtering assum e that F, G, H, K, Q, and R are constant 
and that the steady-state propagation coefficient L is used. Information filter estimation then 
reduces to the following. 


y(Jb|Jb - 1) = Ly(k - l\k - 1) 

y(Ar|/t) = y(k\k — 1) + H T R~ l [z(k) - Ju(k )] 


where 


L = YFY~ l 


and Y = P~ l where P is the solution of the algebraic Riccati equation (3.18). 


(3.26) 

(3.27) 


(3.28) 
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3.4. Comparison of the Kalman and Information Filters 

By comparing the implementation requirements and performance of the Kalman and Infor- 
mation filters, a number of attractive features of the latter are identified: 

• The information estimation equations (3.23) and (3.24) are computationally simpler than 
the state estimation equations (3.10) and (3.11). 

• Although the information prediction equations (3.21) and (3.22) are more complex titan 
equations (3.8) and (3.9), prediction depends on a propagation coefficient which is inde- 
pendent of the observations. It is thus again easy to decouple and decentralize. 

« 

• There are no gain or innovation covariance matrices and the maximum dimension of a 
matrix to be inverted is the state dimension. In multisensor systems the state dimension is 
sometimes smaller than the observation dimension, in which case it is preferable to employ 
the Information filter and invert smaller information matrices than use the Kalman filter 

, and invert larger innovation covariance matrices. 

• Initializing the Information filter is much easier than for the Kalman filter. This is 
because information estimates (matrix and state) are easily initialized to zero information. 
However, in order to implement the Information filter a start up procedure is required, 
where the information matrix is set with small non-zero diagonal elements to make it 
invertible. 

4. Preliminary Results 

This section describes some results using a simplified version of the engine estimation (i.e., 
only a few senors and parameters/input biases axe considered). The results clearly demonstrate 
the utility of using the singular values of A as a means of determining the conditioning of the 
corresponding filtering problem and show that the stacking technique can greatly increase the 
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quality of the estimation of the state and input and parameter biases. The models used at each 
operating point consist of 3 states ( xnl , xnh, tmpc) and the sample period used to generate the 
discrete-time models is T = O.OOlsec. In each simulation, the sensor noise was chosen to have 
a very small covariance and no plant noise was simulated. 

4.1. Example 1 

In the first example the operating point is given by (Altitude=0 ft, Mach No.=0, Power 
Code=40). The bias vector ttf, was of dimension 4 and consisted of 2 input biases (za8«, zp2«) 
and the variations of 2 parameters {sedm2, zsw7d). The sensor vector z consisted of 8 mea- 
surements (xn2c, xn25c, t27c, f56c, psl5c, p27c, ps3c, ps56c) . 

Since A € IR llx7 for this configuration, for rank (A) = ncol(A), the first 7 singular values 
must be nonzero. The first 7 singular values are given by 

‘ 6.77e + 2 ' 

9.25e + 1 
6.99e + 0 

1.01e + 0 . (4.1) 

l.OOe + 0 
1.09e - 1 
8.99e - 3 

Notice that not only are these singular values nonzero, the ratio of the minimum singular value 
to the mavimum singular value is 1.33 E - 5. This probably indicates that the filtering problem 
is well conditioned. This conjecture is clearly verified by Figures 1 thru 7 which show that each 
of the state, input biases and parameter variations converge to their actual values. 
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Figure 1: 8 Measurement Estimation - State 1, XNL 



Figure 2: 8 Measurement Estimation - State 2, XNH 
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Figure 6: 8 Measurement Estimation - Efficiency and Flow Scalar 3, ZSW49 
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TtawSapk 

■* 

Figure 7: 8 Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D 
4.2. Examples 2 and 3 

Example 2. As in Example 1, in this example the operating point is again given by (Alti- 
tude^ ft, Mach No.=0, Power Code=40) and the bias vector u k was of dimension 4 and consisted 
of 2 input biases (za8«,zp2fc) and the variations of 2 parameters (sedml, zsw7d). However, 
now the sensor vector z consisted of only 6 measurements (xn2c,xn25c,t27c,t56c,p27c,ps56c). 
(Example 1 measurements consisted of these 6 measurements plus psl5c and ps3c.) 

Since A € IR 9x7 for this configuration, for rank(A) = ncol(A), the first 7 singular values 

must be nonzero. The first 7 singular values are given by 

' 6.73e + 2 ' 

9.19e + 1 
3.22e + 0 

l.OOe + 0 . (4.2) 

l.OOe + 0 
7.62e - 3 
3.06e - 5 

L J 

Notice that each of these singular values is nonzero, but the ratio of the minimum singular 
value to the maximum singular value is 4.551? — 8 which a factor of 3.41? — 3 smaller than it 
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whs for Example 1. This wider spread may indicate that the filtering problem is ill-conditioned. 
The "unstacked” curves in Figures 8 thru 14 clearly indicate that this is the case since several 
of the states, input biases and parameter variations fail to converge to their actual values. 


Example 3. An attempted remedy to this filter failure is to use data from at least one 
additional operating point. Let us consider one additional operating point described by (Alti- 
tmle=20. 000ft, Mach No.=0.5,Power Code=40) with the same 6 sensors and 4 input biases and 
parameter variations as above. A corresponding to the 2 operating points is 18 x 10. Hence, for 
rnvk(A) = ncol{A). the first 10 singular values must be nonzero. The first 10 singular values 


are given by 

‘ 9.78e + 2 
1.32c + 2 
3.55e + 1 
l.OOe + 0 
l.OOe + 0 
l.OOe + 0 
l.OOe + 0 
1.32e - 1 
4.96e - 3 
1.95e - 3 


(4.3) 


Notice that each of these singular values is nonzero and the ratio of the minimum singular 


value to the maximum singular value is 1.99E — 6 which is a factor of 1.5E-1 smaller than in 
Example 1 but a factor of 4.39 E 4- 1 larger than in Example 2. Since Example 1 converged to 
the correct states, input biases and parameter variations, it is logical to infer that this spread 
indicates that this filtering problem is well conditioned. The “stacked” curves in Figures 8 thru 
14 clearly indicate that this is the case. With the exception of state 3 (see Figure 10) each of 


the states, input biases and parameter variations converged to their actual values. 


NASA/CR— 2006-2 14244 


135 



0 200 400 600 800 1000 1200 1400 1600 1800 200C 

TVn»tapk 


Fieure 8: 6 Measurement Estimation - State 1. XNL 



Figure 9: 6 Measurement Estimation - State 2, XNH 
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Figure 10: 6 Measurement Estimation - State 3, TMPC 



Figure 11: 6 Measurement Estimation - Sensed Input 1, ZA8 
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Figure 12: 6 Measurement Estimation - Environmental Bias 2, ZP2 



Figure 13: 6 Measurement Estimation - Efficiency and Flow Scalar 3, ZSW49 
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Figure 14: 6 Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D 

5. Larger Scale Simulation 

In this section we consider a parameter set more realistic than that considered in the previous 
section. As detailed in the tables of the next subsection, it was assumed that 8 sensors were 
available to estimate 3 states and 15 parameters. The Kalman and Information filters were not 
able to estimate the remaining 8 parameters no matter how many operating points were used. 
Hence, these 8 parameters are not considered in. this simulation. 

Linear models were constructed for 59 operating points in the flight envelope. These oper- 
ating points are described in Appendix 1. A sampling period of T=0.2 sec was used for the 
discretization of these linear models. In each simulation the actual values of the parameter set 
were calculated without white process noise, w(k). The white sensor noise, v(Ar), was simulated 
using the MATLAB function RAND to sample a uniform distribution between 0 and 1. In the 
Kalman filter algorithm the covariance, Q(k), of the process noise and the covariance, R{k), 
of the sensor noise were chosen to be diagonal. The diagonal elements of of R(k) were chosen 
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to be equal to 0.1. The diagonal elements of Q(k) were chosen such that the three diagonal 
elements corresponding to the three plant states were 0.1, the variances of the 3 directly sensed 
inputs, ZWF36, ZA8, ZA16 (parameters 1,2, and 3 in Table 4) were 0.001, and the remaining 
diagonal elements were 0.0001. 

Since 15 parameters were used in the simulation, the dimensions of A, defined by (2.9) 
differed from those tabulated in Table 1. The dimensions for combinations of one, two, and 
three operating points are given in Table 2. For combinations of two and three operating points 
the inverse of the condition number of A, cond(A), defined by 

amd(A) = ^£>1, 

°mtn 

was used to determine how well the matrix satisfied the full column rank condition. In particular 
the larger the value of amd(A)~ l the better a combination of two or three operating points. 
(Recognize that cond(.4.) -1 < 1.) All combinations of two and three operating points were 
created from the 59 available operating points and the corresponding cond{A)~ l were computed. 


nr 

nrow(.A) 

ncol(.A) 

i 

11 

18 

2 

22 

21 

3 

33 

24 


Table 2: Number of Rows and Columns of A 
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The largest 100 cond(A)~ l for combinations of two and three operating points are listed in 
Appendix 1. Comparison of the largest amd(A)~ l for combinations of two operating points, 
6.68E-Q7, to that of the largest cond(A)~ l using three operating points, 2.72E-06, reveals tnat 
the latter number is larger by a factor of 4. As seen in Appendix 1, the best combination of 
two operating points is (39, 40), while the best combination of three operating points is (7, 51, 
58). Of the total 59 operating points considered only 43 were used in the top 100 combinations 
of three, while all were used in the top 100 combinations of two. 

Except in the cases where the Kalman filter diverged, the estimates of the Kalman and 
Information filters was quite close. Hence, in the tabular results reported in this section (see 
Tables 5, 6 and 7), we report only the statistics associated with parameter estimation using 
the Information filter. The statistics for both the Kalman and Information filters are given in 
Tables 12 thru 18 of Appendices 2 thru 8. 

5.1. Sensors and Parameters 


Sensor Number 

Sensor Name 

1 

XN2C - Low pressure rotor speed 

2 

XN25C - High pressure rotor speed 

3 

T27C - Core driven fan stage hub temperature 

4 

T56C - LP Turbine exit temperature 

5 

PS15C - Bypass duct static pressure 

6 

P27C - Core driven fan stage hub total pressure 

7 

PS3C - Compressor discharge static pressure 

8 

PS56C - LP Turbine frame static pressure 


Table 3: Sensors Available for Simulation 
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Parameter Number 

Parameter Name 

Parameter Value 

1 

ZWF36 - Combustor fuel flow 


2 

ZA8 - Nozzle area 

0.00E+00 

3 

ZA16 - Rear VABI area 

O.OOE+OO 

4 

ZSW2 - Fan flow scalar 

3.00E-03 

5 

SEDM2 - Fan efficiency scalar 

3.00E-03 

6 

ZSW7D - Compressor tip flow scalar 

3.00E-03 

7 

SEDM7D - Compressor tip efficiency scalar 

3.00E-03 

8 

ZSW27 - Compressor hub flow scalar 

3.00E-03 

9 

SEDM27 - Compressor hub efficiency scalar 

3.00E-03 

10 

ZSW41 - HP Turbine flow scalar 

3.00E-03 

11 

ZSE41 - HP Turbine efficiency scalar 

3.00E-03 

12 

ZSW49 - LP Turbine flow scalar 

3.00E-03 

13 

ZSE49 - LP Turbine efficiency scalar 

o not? no « 

U.UUU'VU 1 

14 

ZWB3 - Compressor discharge bleed flow 

0.00E+00 

15 

ZPWXH - Customer power extraction (core rotor) 

0.00E+00 


Table 4: Parameter Set to Be Estimated in Simulation 


5.2. Simulation Results 

Below several simulation and tabular results for nine cases are given. Complete sets of the 
simulation and tabular results are presented for the first seven cases in Appendices 2 thru 8. 
For the first four cases each of the parameters was initialized to zero in the filtering. These four 
cases respectively illustrate estimation using 1) the best combination of three operating points, 
2) an arbitrary combination of three operating points, 3) the best combination of two operating 
points, and 4) a single operating point from the best combination of three. In each of these 
cases the parameter sets converge nicely, although some results do have better convergence 
than others. 

In the next three cases each of the parameters was initialized to 0.5. Estimation results are 
given for the cases described above, with the exception of the arbitrary combination of three 
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operating points. When the parameters are initialized to 0.5, it is seen that the parameters 
do not converge to their true values using the best combination of two operating points. The 
results are even worse using the single operating point. 

The final two cases illustrate that the Kalman filter sometimes diverged while the Informa- 
tion filter converged. Hence, the Information filter is seen to be more numerically robust. 

5.2.1. Case 1: Estimation Using The Best Combination of Three • Parameters 
Initialized to Zero 

The first case considered was the combination of three operating points (7,51.58) with 
the largest cond{A)~ l . Table 5 and Figures 15 thru 17 show that the data from these three 
operating points are sufficient to estimate the parameter set successfully. A complete set of 
figures illustrating the estimation is given in Appendix 2. 



Figure 15: Measurement Estimation - Sensed Input 1, ZWF36 (Case 1) 
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Figure 16: Measurement Estimation - Sensed Input 2, ZA8 (Case 1) 



Figure 17: Measurement Estimation - Sensed Input 3, ZA16 (Case 1) 
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5.2.2. Case 2: Estimation Using Arbitrary Combination of Three - Parameters 
Initialized to Zero 

An arbitrary combination of three operating points was chosen (1,30,57) to be used to 
estimate the parameter set. For this case cond(A )~ 1 = 8.73E-7, which is a factor of 3 smaller 
than the previous case. Table 5 and Figures 18 thru 20 show that, as expected, the estimation 
process is not as accurate as for the best combination of three operating points. However, the 
results are still reasonable. A complete set of figures illustrating the estimation is given in 
Apppndix 3. 


i Hm tor MM 1 Mtof towa < 



Figure 18: Measurement Estimation - Sensed Input 1, ZWF36 (Case 2) 
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Case 1 


Case 2 


Parameter 

Value 

Mean 

Variance 

Mean 

Variance 

ZWF36 

3.00E-03 

3.24E-02 

5.20E-03 

3.24E-02 

5.20E-03 

ZA8 

0.00E+00 

5.20E-03 

6.20E-04 

-4.75E-04 

6.20E-04 

ZA16 

0.00E+00 

2.48E-02 

2.97E-03 

2.48E-02 

2.97E-03 

ZSW2 

3.00E-03 

2.99E-03 

2.06E-07 

2.99E-03 

2.06E-07 

SEDM2 

3.00E-03 

2.99E-03 

1.06E-07 

2.99E>03 

1.06E-07 

ZSW7D 

3.00E-03 

2.97E-03 

6.46E-07 



SEDM7D 

3.00E-03 

2.95E-03 

1.29E-06 

2.95E-03 

1.29E-06 

ZSW27 

3.00E-03 

3.02E-03 

3.89E-08 

3.02E-03 

3.89E-08 

SEDM27 

3.00E-03 

3.01E-03 

7.10E-08 

3.01E-03 

7.10E-08 

ZSW41 

3.00E-03 

3.01E-03 

5.05E-08 

3.01E-03 

5.05E-08 

ZSE41 

3.00E-03 

2.99E-03 

1.62E-07 

2.99E-03 

1.62E-07 

ZSW49 

3.00E-03 

3.02E-03 

2.12E-07 

3.02E-03 

2.12E-07 

ZSE49 

3.00E-03 

3.02E-03 

2.00E-07 

3.02E-03 

2.00E-07 

ZWB3 

0.00E+00 

1.56E-04 

2.11E-05 

1.56E-04 

2.11E-05 

ZPWXH 

O.OOE+OO 

-2.21E-02 

6.81E-04 

-2.21E-02 

6.81E-04 


Table 5: Means and Variances of Parameters Estimated Using the Information Filter - Cases 1 
and 2 


5.2.3. Case 3: Estimation Using Best Combination of Two - Parameters Initialized 
to Zero 


Next consider the combination of two operating points (39,40) with the largest cond{A)~ l 
(6.68E-07) which is less than either of the combinations of three considered previously. Table 6 
and Figures 21 thru 23 reveal that, as expected, the parameter estimation has worse convergence 
properties than the previous cases which used three operating points. A complete set of figures 
illustrating the estimation is given in Appendix 4. 
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Figure 21: Measurement Estimation - Sensed Input 1, ZWF36 (Case 3) 



Figure 22: Measurement Estimation - Sensed Input 2, ZA8 (Case 3) 
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Figure 23 : Measurement Estimation - Sensed Input 3, ZA16 (Case 3) 

5 . 2 . 4 . Case 4 : Estimation Using Single Operating Point - Parameters Initialized 
to Zero 

It was determined by simulation that the single operating point from the best combination 
of three (7,51,58) that gives the best estimation results is 58. Since for one operating point 
A € IR llxI8 such that nrow(A) < ncol(A), it is possible for the estimation process to converge 
to the wrong parameter values. However, Table 6 and Figures 24 thru 26 illustrate that the 
estimation process with a single operating point can yield good results. This is because for this 
and all previous cases the parameter set was initialized to zero for the filtering, such that the 
initial parameter estimates were equal or close to their actual values. It will be shown below 
that the estimation process for a single operating point yields poor estimates with different 
initial parameters. A complete set of figures illustrating the estimation is given in Appendix 5. 
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Figure 24: Measurement Estimation - Sensed Input 1, ZWF36 (Case 4) 



Figure 25: Measurement Estimation - Sensed Input 2, ZA8 (Case 4) 
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Figure 26: Measurement Estimation - Sensed Input 3, ZA16 (Case 4) 




Case 3 


Case 4 


Parameter 

Value 

Mean 

Variance 

Mean 

Variance 

ZWF36 

3.00E-03 

6.57E-03 

1.52E-02 

7.26E-05 

2.30E-08 

ZA8 

0.00E+00 

3.90E-02 

1.70E-02 

-2.70E-03 

2.00E-05 

ZA16 

0.00E+00 

1.12E-01 

9.10E-03 

1.23E-03 

6.11E-06 

ZSW2 

3.00E-03 

2.94E-03 

1.56E-05 

2.95E-03 

2.98E-06 

SEDM2 

3.00E-03 

2.90E-03 

6.52E-06 

2.94E-03 

4.05E-06 

ZSW7D 

3.00E-03 

3.21E-03 

4.56E-05 

2.94E-03 

6.33E-06 

SEDM7D 

3.00E-03 

9.21E-04 

9.42E-05 

4.03E-03 

5.05E-05 

ZSW27 

3.00E-03 

2.92E-03 

1.95E-06 

3.15E-03 

4.97E-06 

SEDM27 

3.00E-03 

3.11E-03 

2.38E-06 

3.45E-03 

4.67E-07 

ZSW41 

3.00E-03 

2.98E-03 

1.29E-06 

3.05E-03 

5.77E-07 

ZSE41 

3.00E-03 

3.05E-03 

1.29E-06 

3.42E-03 

1.94E-06 

ZSW49 

3.00E-03 

3.01E-03 

7.03E-06 

1.83E-03 

8.78E-07 

ZSE49 

3.00E-03 

3.05E-03 

1.48E-06 

4.64E-04 

5.59E-07 

ZWB3 

O.OOE-fOO 

KItEElgl 

5.18E-04 


6.78E-07 

ZPWXH 

0.00E+00 

1.08E-02 

1.76E-04 

-1.47E-06 

6.14E-12 


Table 6: Means and Variances of Parameters Estimated Using the Information Filter - Cases 3 
and 4 
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5.2.5. Case 5: Estimation Using Best Combination of Three - Parameters Initial- 
ized to 0.5 

This case parallels that described in subsection 5.2.2 for the combination (7,51,58) with the 
exception that each of the parameters was initialized to 0.5. Table 7 and Figures 27 thru 29 
show that the parameters converge to their actual values despite the poor initialization. A 
complete set of figures illustrating the estimation is given in Appendix 6. 


Figure 27: Measurement Estimation - Sensed Input 1, ZWF36 (Case 5) 
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5.2.6. Case 6: Estimation Using Best Combination of Two - Parameters Initialized 
to 0.5 

This case parallels that described in subsection 5.2.3 for the combination (39,44) with the 
exception that each of the parameters was initialized to 0.5. Table 7 and Figures 30 thru 32 
illustrate that the parameters converge slowly and do not always converge to their true values. 
A complete set of figures illustrating the estimation is given in Appendix 7. 
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Figure 30: Measurement Estimation - Sensed Input 1, ZWF36 (Case 6) 
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Figure 31: Measurement Estimation - Sensed Input 2, ZA8 (Case 6) 
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Figure 32: Measurement Estimation - Sensed Input 3, ZA16 (Case 6) 
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Case 5 


Case 6 


Parameter 

Value 

Mean 

Variance 

Mean 

Variance 

ZWF36 

3.00E-03 

4.34E-02 

1.19E-02 

1.15E-01 

1.65E-02 

ZA8 

0.00E+00 

-1.94E-02 

7.97E-04 

1.53E-01 

MESraStKlI 

ZA16 

0.00E+00 

1.19E-02 

2.51E-03 

2.74E-01 

2.25E-02 

ZSW2 

3.00E-03 

2.95E-03 

2.07E-07 

3.93E-03 

3.04E-05 

SEDM2 

3.0QE-03 

2.99E-03 

1.07E-07 

3.30E-03 

9.50E-06 

ZSW7D 

3.00E-03 

2.96E-03 

6.49E-07 

1.77E-03 

5.16E-05 

SEDM7D 

3.00E-03 

2.81E-03 

1.33E-06 



ZSW27 

3.00E-03 

3.01E-03 

3.85E-08 

3.39E-03 

6.05E-06 

SEDM27 

3.00E-03 

2.97E-03 

7.22E-08 

2.91E-03 

2.60E-06 

ZSW41 

3.00E-03 

3.00E-03 

4.90E-08 

3.15E-03 

1.51E-06 

ZSE41 

3.00E-03 

3.10E-03 

1.68E-07 

3.06E-03 

3.64E-06 

ZSW49 

3.00E-03 

3.00E-03 

2.11E-07 

2.82E-03 

1.26E-05 

ZSE49 

3.00E-03 

3.01E-03 

1.93E-07 

3.12E-03 

1.93E-06 

ZWB3 

O.OOE+OO 

1.49E-04 

2.33E-05 

3.17E-03 

1.65E-03 

ZPWXH 

O.OOE+OO 

1.72E-01 

2.56E-02 

4.67E-01 

1.75E-03 


Table 7: Means and Variances of Parameters Estimated Using the Information Filter - Cases 5 
and 6 


5.2.7. Case 7: Estimation Using Single Operating Point from Best Combination 
of Three - Parameters Initialized to 0.5 


This case parallels that described in subsection 5.2.4 with the exception that each of the 
parameters was initialized to 0.5. Table 7 and Figures 33 thru 35 show that the parameters do 
not converge to their true values. A complete set of figures illustrating the estimation is given 
in Appendix 8. 
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Figure 33: Measurement Estimation - Sensed' Input 1, ZWF36 (Case 7) 
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Figure 34: Measurement Estimation - Sensed Input 2, ZAo (Case 7) 
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Case 7 

■■■ 

Parameter 

Value 

Mean 

Variance 

ZWF36 

3.00E-03 

4.94E-01 

2.13E-08 

ZA8 

0.00E+00 

7.66E-01 

1.86E-05 

ZA16 

O.OOE+OO 

3.98E-01 

5.70E-06 

ZSW2 

3.00E-03 

2.09E-02 

3.26E-06 

SEDM2 

3.00E-03 

4.03E-03 

4.61E-06 

ZSW7D 

3.00E-03 

7.02E-03 

6.83E-06 

SEDM7D 

3.00E-03 

1.14E-01 

4.68E-05 

ZSW27 

3.00E-03 

4.42E-02 

5.46E-06 

SEDM27 

3.00E-03 

-1.55E-01 

5.27E-07 

ZSW41 

3.00E-03 

1.32E-02 

5.64E-07 

ZSE41 

3.00E-03 

-3.99E-02 

1.91E-06 

ZSW49 

3.00E-03 

2.22E-01 

8.95E-07 

ZSE49 

3.00E-03 

4.22E-01 

5.66E-07 

ZWB3 

0.00E+00 

5.61E-01 

6.27E-07 

iMsmmmism 

5.00E-01 

5.67E-12 


Table 8: Means and Variances of Parameters Estimated Using the Information Filter - < 
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5.2.8. Estimation Illustrating Kalman Filter Divergence and Information Filter 
Convergence 

One of the advantages of the Information filter over the Kalman filter (as discussed in 
Section 3.4) is the fact that there are no gain W(k) nor innovation covariance S(k) matrices 
involved in its algorithm, whereas these two matrices are central to the Kalman filter as shown in 
(3.10) - (3.12). From (3.11) it is clear that the gain W(k) depends on the predicted covariance 
P(fcj& — 1), observation matrix H(k) and the innovation covariance S(k.) The innovation 
covariance in turns depends on H(k), P(k\k - 1) and R(k). This means the critical variable 
in both W(k) and S(k) is the predicted covariance P(fc|A: - 1). This predicted covariance is 
computed from the previous covariance P(k — l|fc — 1) given by (3.9). This means that any 
problems associated with an ill-conditioned covariance P(k\k), such as singularity or divergence, 
will affect the Kalman filter. 

On the other the hand, the Information filter does not involve W(A) nor S(k) at all. Al- 
though the definition of Fisher information given in (3.19) defines the Information matrix Y(k\k) 
as the inverse of the covariance P(k\k), a careful inspection of the Information filter algorithm 
shows that it is not necessary to continuously compute and use P(k\k). Consequently any 
problems associated with an ill-conditioned P(k\k), S(k) or W(k) won’t affect the Information 
filter. 

In particular, the concerns about the singularity of P(k\k) which are raised in the literature 
[12], are of no consequence to the Information filtering techniques discussed in this report. The 
only covariance required is the initial one, i.e P(0|0), to obtain y(0|0). The Information filter 
algorithm avoids the singularity of P(0|0) by using a start-up procedure for initialization, where 
the initial covariance is equated to a product of a constant a (a large positive number such as 
100 or 1000) and the the process noise covariance Q(k) as follows: 

P(0|0) = oQ(0) (5.1) 
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r(o|o) = P(ojo )- 1 


(5.2) 


The idea is to have a start-up procedure where the initial information matrix is set with 
small non-zero diagonal elements [7]. This is in keeping with the fact that initializing the 
Information filter is much easier than initializing the Kalman filter, because the Information 
estimates (matrix and state) are easily initialized to zero information. From (3.21) - (3.25) it 
is clear that once the Information filter is started it does not require nor use P(Jfc|&) and hence 
the singularity of P(ArjA) is of no relevance. 

Another form of an ill-conditioned P(fc|A:) is when the covariance matrix diverges, that 
is, its elements grow without bound towards infinity. It is seen from (3.13) that when P(k\jc) 
diverges, (due to numerical ill-conditioning) the innovation covariance matrix S(k) also diverges. 
Consequently, from (3.10) the state estimate m(Ar|Ar) similarly diverges, producing the divergence 
of the Kalman filter algorithm. On the other hand, since the terms P{k\k), W(k ) and S(k) do 
not appear in the Information filter algorithm, it does not suffer from this divergence problem 
and hence successfully converges. 

Figures 36 thru 41 show two cases in which the Kalman filter diverges due to the divergence 
of P{k\k) and S(k), while the Information filter estimates successfully. Figures 36 thru 38 
correspond to the combination (10,46,58) with cond(A)~ l = 1.45E-7, while Figures 39 thru 
41 correspond to the combination (1,2,3) with cond(A) -1 = 9.43E-8. These results effectively 
illustrate the advantages of the Information filter over the Kalman filter in situations where 
combinations of operating points lead to ill-conditioned covariance matrices. 
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Figure 36: Measurement Estimation - Sensed Input 1, ZWF36 (Case 8) 



Figure 37: Measurement Estimation - Sensed Input 2, ZA8 (Case 8) 
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Figure 38: Measurement Estimation - Sensed Input 3, ZA16 (Case 8) 
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Figure 39: Measurement Estimation - Sensed Input 1, ZWF36 (Case 9) 
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Figure 40: Measurement Estimation - Sensed Input 2, ZA8 (Case 0) 
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Figure 41: Measurement Estimation - Sensed Input 3, ZA16 (Case 9) 
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6. Conclusions 


This report has studied the limitations of estimating the state and parameter and input 

biases at one operating point in the flight envelope. It was seen that if the number of sensors 

is less than the number of biases to be estimated, then it is essential to use data from more 

than one operating point to perform the estimation. A procedure for calculating the minimum 

number of operating points and for determining an appropriate combination of operating points 

was given in terms of a rank test of a matrix A that is a function of the system matrices of 

each of the operating points. In addition, a “stacked” model structure was proposed as the 

basis of the estimation. Practical implementation of the estimation would then require storing 

« • 

data from several operating points and some preliminary results of that implementation have 
been given. The results in this report validate the use of the matrix A in predicting estimation 
performance. 

Given a realistic set of 59 operating points from the flight envelope, all combinations of three 
and two were computed as well as the condition numbers of A for each case, which determines 
essentially how well A meets the rank condition. It was shown that this number has a strong 
bearing on how well the parameter set is estimated. The data from three appropriately chosen 
operating points was shown both theoretically and via simulation to be adequate for estimating 
t he parameter set. 

The results reported here clearly indicate that with the “stacked” configuration the use of 
Information filter will be more numerically robust than Kalman filter. While the Kalman filter 
sometimes diverged in our simulation, the Information filter always converged. 

Finally, it should be noted that although the results developed here assume that exclu- 
sively steady-state data is available for estimation. In practice, a combination of transient and 
steady-state data will be available. When transient data is used, the analysis involving A is 
not strictly valid as presented here. This is because the problem becomes significantly more 

59 
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complex. In particular, the plant matrices, A , B, C, etc..., vary with time and the steady-state 
assumptions do not hold. In this case, standard Kalman and Information filtering are also not 
adequate. Instead, it is necessary to use Extended Kalman and Information filtering. However, 
the “stacked” configuration should still be a useful paradigm for obtaining better estimates. 
Each level of the stack will now correspond to a particular set of transient and steady-state 
data where each of these data sets is of the same length. Unlike the assumptions in this report, 
the data sets will in general not represent data about one operating point. Extended Kalman 
or Information filtering must then be used to perform the estimation. 

The question remains as to which data sets will be adequate? Although the analysis with A 

* 

is not strictly valid since it assumes steady-state data, it may still be useful. Using engineering 
judgment, it is possible to associate with each set of (transient and steady-state) data a finite 
number of plant matrix sets, representing the system at various points in the tiignt regimes 
represented by the data set. Hence, we can associate with each combination of data sets, finite 
sets of A matrices. It is conjectured that evaluation of these sets of A matrices, will provide 
valuable insight into which combinations of data sets will yield the most accurate estimates. 
However, the details of the exact procedures to use involving transient data must be developed 
and could be the subject of further research. 
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Appendix 1: Tables of Operating Points and Combination 
Rankings 


The first table lists each of the 59 operating points along with the corresponding altitude, 
mach number, ambient temperature and power code. The second table lists the top 100 com- 
binations of two operating points, ranked according to cond{A)~ 1 . The third table lists the top 
100 combinations of three opertaing points, ranked according to cond(A)~ l . 


Operating Point 

Altitude 

Mach No 

Amb. Temperature 

Power Code 

1 

0 

0 

h 0 

50 

2 

0 

0 

0 

47.5 

3 

0 

0 

0 

45 

4 

0 

0 

0 

42.5 

5 

0 

0 

0 

40 

lam 

0 

0 

0 

35 

mmmmm 

0 

0 

0 

27.5 

8 


0 

0 

25 

9 

0 

0 

0 

22.5 

10 

0 

0 

0 


11 

40000 

0.4 

0 

50 


40000 

0.4 

0 

47.5 


40000 

0.4 

0 

45 i 

14 

40000 

0.4 

0 

42.5 


40000 

0.4 

0 

40 

16 

40000 

0.4 

0 

37.5 


40000 

0.4 | 

0 

32.5 

18 

40000 

0.4 

0 

30 

19 

40000 

0.4 

0 

27.5 


40000 

0.4 

0 

22.5 

21 

40000 

0.4 

0 

21 

22 

0 

0 

-100 

50 

23 

0 

0 

-50 

50 

24 

45000 

1.725 

0 


25 

50000 

1.725 

0 

50 
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Table 9: Operating Points Used From The Flight Envelope 
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Combination 

Operating Points 

Cond(A)~ l 

1 

39 44 

6.6800E-07 

2 

18 47 

3.0446E-07 

3 

39 51 

2.8830E-07 

4 

39 47 

2.7890E-07 

5 

7 39 

2.5803E-07 

6 

13 50 

2.5368 E-07 

7 

39 57 

2.3908E-07 

8 

34 40 

2.2053E-07 

9 

33 39 

2.2027E-07 

10 

37 49 

2.1896E-07 

11 

4 39 

2.1397E-07 

12 

24 37 

2.1244E-07 

13 

33 37 

2.1188E-07 

14 

21 39 

2.0969E-07 

15 

13 44 

2.0730E-07 

16 

30 39 

2.0536E-07 

17 

40 47 

2.0015E-07 

18 

13 18 

1.9843E-07 

19 

21 40 

1.9829E-07 

20 

31 39 

1.9634E-07 

21 

13 41 

1.8841E-07 

22 

39 46 

1.8527E-07 

23 

19 47 

1.8162E-07 

24 

18 34 

1.7726E-07 

25 

3 39 

1.7034E-07 

26 

38 47 

1.6883E-07 

27 

37 43 

1.6695E-07 

28 

27 39 

1.6140E-07 

29 

28 55 

1.5912E-07 

30 

6 51 

1.5757E-07 

31 

48 57 

1.5683E-07 

32 

21 56 

1.5595E-07 

33 

13 21 

1.5420E-07 

34 

13 59 

1.5152E-07 

35 

25 37 

1.5131E-07 
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| Combination 

Operating Points 

Cond(A)~ l 

i 36 

26 39 

1.5072E-07 

37 

21 27 

1.4873E-07 

38 

37 46 

1.4829E-07 

! 39 

39 43 

1.4638E-07 

j 40 

8 39 

1.4081E-07 

1 41 

32 39 

1.4035E-07 

42 

39 54 

1.3693E-07 

j 43 

25 39 

1.3606E-07 

44 

12 39 

1.3408E-07 

| 45 

24 39 

1.3362E-07 

46 

2 39 

1.3310E-07 

| 47 

11 39 

1.3291E-07 

48 

39 59 

1.3207E-07 

49 

1 39 

1.3111E-07 

50 

29 39 

1.2847E-07 

51 

39 48 

1.2755E-07 

52 

13 39 

1.2709E-07 

53 

35 44 

1.2691E-07 

54 

39 49 

1.2670E-07 

55 

15 55 

1.2604E-07 

56 

18 53 

1.2324E-07 

57 

7 55 

1.2080E-07 

58 

30 37 

1.1927E-07 

59 

13 25 

1.1831E-07 

60 

20 47 

1.1783E-07 

61 

14 39 

1.1715E-07 

62 

13 24 

1.1662E-07 

63 

35 39 

1.1626E-07 

64 

7 19 

1.1539E-07 

65 

18 58 

1.1420E-07 

66 

36 55 

1.1346E-07 

67 

34 38 

1.1330E-07 

68 

32 37 

1.1296E-07 

69 

16 55 

1.1130E-07 

70 

39 42 

1.1129E-07 
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6 10 

9.7203E-08 

38 51 

9.6950E-08 

18 42 

9.6464E-08 

39 58 

9.4706E-08 

12 37 

9.4689E-08 

17 51 

9.4608E-08 


Table 10: Laregest 100 Singular Values and Corresponding Combinations of Two Operating 
Points 
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Combination 

Operating Points 

WEESBHMk 

1 

7 51 58 

2.7247E-06 

2 

18 44 47 

2.5709E-06 

3 

8 42 56 

2.4606E-06 

4 

28 31 42 

2.4603E-06 

5 

18 31 44 

2.3981E-06 

6 

18 44 51 

2.3901E-06 1 

7 

9 13 18 

2.3594E-06 

8 

8 18 47 

2.3505E-06 

9 

13 47 51 

2.3488E-06 

10 

19 44 47 

2.3435E-06 

11 

28 34 42 

2.3365E-06 

12 

12 47 51 

2.3321E-06 

13 

18 34 38 

2.3309E-06 

14 

8 18 31 

2.3271E-06 

15 

18 38 47 

2.3223E-06 

16 

27 31 48 

2.3223E-06 

17 

17 34 44 

2.3213E-06 

18 

27 47 51 

2.3137E-06 

19 

17 44 47 

2.2952E-06 

20 

24 51 55 

2.2924E-06 

21 

18 27 51 

2.2833E-06 

22 

8 18 34 

2.2825E-06 

23 

8 48 56 

2.2803E-06 

24 

10 18 19 

2.2799E-06 

25 

9 12 18 

2.2755E-06 

26 

27 48 51 

2.2672E-06 

27 

17 31 44 

2.2557E-06 

28 

9 18 43 

2.2544E-06 

29 

42 47 51 

2.2536E-06 

30 

37 44 51 

2.2509E-06 

31 


2.2498E-06 

32 

9 14 18 

2.2447E-06 

33 

18 21 47 

2.2338E-06 

34 

11 47 51 

2.2319E-06 

35 

9 18 27 

2.2230E-06 
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Combination 

Operating Points 

Cond(A) 1 

36 

13 18 57 

2.2194E-06 

37 

27 31 58 

2.2168E-06 

38 

9 16 18 

2.2041E-06 

39 

18 20 51 

2.1990E-06 

40 

18 21 34 

2.1975E-06 

41 

28 31 58 

2.1954E-06 

42 

18 38 51 

2.1945E-06 

43 

7 47 59 

2.1933E-06 

44 

26 47 51 

2.1873E-06 

45 

10 18 44 

2.1829E-06 

46 

17 20 47 

2.1824E-06 

47 

19 44 51 

2.1822E-06 

48 

7 47 49 

2.1821E-06 

49 

24 28 31 

2.1793E-06 

50 

17 20 34 

2.1779E-06 

51 

25 28 31 

2.1746E-06 

52 

13 32 51 

2.1728E-06 

53 

27 47 59 

2.1701E-06 

54 

8 56 58 

2.1681E-06 

55 

9 15 18 

2.1668E-06 

56 

34 37 44 

2.1643E-06 

57 

27 51 53 

2.1629E-06 

58 

10 19 20 

2.1628E-06 

59 

8 10 18 

2.1624E-06 

60 

47 56 59 

2.1619E-06 

61 

11 14 30 

2.1612E-06 

62 

9 44 47 

2.1562E-06 

63 

28 31 48 

2.1559E-06 

64 

35 38 44 

2.1542E-06 

65 

11 15 30 

2.1522E-06 

66 

8 13 18 

2.1515E-06 

67 

7 18 31 

2.1502E-06 

68 

19 27 31 

2.1456E-06 

69 

8 15 43 

2.1454E-06 

70 

17 44 51 

2.1448E-06 
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Table 11: Largest 100 Singular Values and Corresponding Combinations of Three Operating 
Points 
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Appendix 2: Complete Data for Estimation Using the 
Best Combination of Three - Parameters Initialized to 
Zero (Case 1) 

Below presents 1) a complete set of figures illustrating the estimation process. 2) a ta- 
ble describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters, and 3) a complete list of the singular values of A. 


Katmsn and Worm — tor tap u* i utang tar** op r— n q puns 



Figure 42: Measurement Estimation - Sensed Input 1, ZWF36 (Case 1) 
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Figure 43: Measurement Estimation - Sensed Input 2, ZA8 (Case 1) 
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Figure 44: Measurement Estimation - Sensed Input 3, ZA16 (Case 1) 
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Figure 45: Measurement Estimation - Efficiency and Flow Scalar 1, ZSW2 (Case 1) 
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Figure 47: Measurement Estimation- Efficiency and Flow Scalar 3, ZSW7D (Case 1) 
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Figure 50: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case 1) 
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Figure 51: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 1) 
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Figure 52: Measurement Estimation - Efficiency and Flow Scalar 8, ZSE41 (Case 1) 
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Figure 55: Measurement Estimation - Bleed 1, ZWB3 (Case 1) 
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Figure 56: Measurement Estimation - Power Extraction 1, ZPWXH (Case 1) 
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Parameter 

Actual Value 

Mean:KF 

Mean:IF 

Variancce:KF 

Variance:IF 

ZA8 

0.00E+00 

-1.74E-03 

-4.75E-04 

9.56E-04 


ZA16 

0.00E+00 

2.41E-02 

2.48E-02 

3.17E-03 


ZSW2 

3.00E-03 

2.99E-03 

2.99E-03 

2.34E-07 


SEDM2 

3.00E-03 

3.00E-03 

2.99E-03 

1.13E-07 

1.06E-07 

ZSW7D 

3.00E-03 

2.97E-03 

2.97E-03 

7.38E-07 

6.46E-07 


3.00E-03 

2.95E-03 

2.95E-03 


1.29E-06 

ZSW27 

3.00E-03 

3.02E-03 

3.02E-03 

4.64E-08 


SEDM27 

3.00E-03 

3.00E-03 

3.01E-03 

7.64E-08 


ZSW41 

3.00E-03 

3.01E-03 

3.01E-03 

8.32E-08 

5.05E-08 

ZSE41 

3.00E-03 

3.00EMJ3 

2.99E-03 

1.45E-07 

1.62E-07 

ZSW49 

3.00E-03 

3.01E-03 

3.02E-03 

1.80E-07 

2.12E-07 

ZSE49 

3.00E-03 

3.02E-03 

3.02E-03 

1.73E-07 

2.00E-07 


O.OOE+OO 

2.04E-04 

1.56E-04 

2.11E-05 

2.11E-05 

| ZPWXH 

0.00E+00 

-1.97E-02 

-2.21E-02 

1.21E-03 

6.81E-04 


Table 12: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Cas 
1 ) 
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Here is the list of singular values: 


3.9709588e + 03 
2.8197357e + 03 
1.0810832e + 03 
7.1665676e + 02 
3.9900754e + 02 
3.1146032e + 02 
2.1304143e + 02 
1.6403454c + 02 
1.1763709c + 02 
2.8748288e + 01 
2.3572168e + 00 
1.0517054c + 00 
1.0381604c + 00 
1.0117562c + 00 
9.9682120e - 01 
9.5306144e - 01 
9.4553955e - 01 
3.1604635e - 01 
7.0784597e - 02 
3.2675650e - 02 
2.9271329e - 02 
1.9956495e - 02 
1.4331101e - 02 
1.0819510e - 02 
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Appendix 3: Complete Data for Estimation Using an Ar- 
bitrary Combination of Three - Parameters Initialized to 
Zero (Case 2) 

Below presents 1) a complete set of figures illustrating the estimation process, 2) a ta- 
ble describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters, and 3) a complete list of the singular values of A. 



Figure 57: Measurement Estimation - Sensed Input 1, ZWF36 (Case 2) 
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Figure 58: Measurement Estimation - Sensed Input 2, ZA8 (Case 2) 
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Figure 59: Measurement Estimation - Sensed Input 3, ZA16 (Case 2) 
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Figure 61: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 2) 
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Figure 62: Measurement Estimation- Efficiency and Flow Scalar 3, ZSW7D (Case 2) 



Figure 63: Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D (Case 2) 
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Figure 64: Measurement Estimation- Efficiency and Flow Scalar 5, ZSYV27 (Case 2) 




Figure 65: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case 2 
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Figure 66: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 2) 



Figure 67: Measurement Estimation - Efficiency and Flow Scalar 8, ZSE41 (Case 2) 
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Figure 68: Measurement Estimation - Efficiency and Flow Scalar 9, ZSW49 (Case 2) 
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Figure 69: Measurement Estimation - Efficiency and Flow Scalar 10, ZSE49 (Case 2) 
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Figure 70: Measurement Estimation - Bleed 1, ZWB3 (Case 2) 
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Figure 71: Measurement Estimation - Power Extraction 1, ZPWXH (Case 2) 
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Parameter 

Actual Value 

Mean:KF 

Mean:IF 

Variancce:KF 

Variance:IF 

ZA8 

O.OOE+OO 

-1.74E-03 

-4.75E-04 

9.56E-04 

6.20E-04 

ZA16 

0.00E+00 

2.41E-02 

2.48E-02 

3.17E-03 

2.97E-03 

IbEbSSkH 

3.00E-03 

2.99E-03 

2.99E-03 

2.34E-07 

2.06E-07 

SEDM2 

3.00E-03 

3.00E-03 

2.99E-03 

1.13E-07 

1.06E-07 

ZSW7D 

3.00E-03 

wmmmwmmm 

7.38E-07 

6.46E-07 

SEDM7D 

3.00E-03 


9.48E-07 

1.29E-06 

ZSW27 

3.00E-03 

3.02E-03 

3.02E-03 

4.64E-08 

3.89E-08 

SEDM27 

3.00E-03 

3.00E-03 

3.01E-03 

7.64E-08 

7.10E-08 

ZSW41 

3.00E-03 

3.01E-03 

3.01E-03 

8.32E-08 

5.05E-08 

ZSE41 

3.00E-03 

3.00E-03 

2.99E-03 

1.45E-07 

1.62E-07 

ZSW49 

3.00E-03 

3.01E-03 

3.02E-03 

1.80E-07 

2.12E-07 

ZSE49 

3.00E-03 

3.02E-03 

3.02E-03 

1.73E-07 

2.00E-07 

ZWB3 

0.00E+00 

2.04E-04 

1.56E-04 

2.11E-05 

2.11E-05 

ZPWXH 

0.00E+00 

-1.97E-02 

-2.21E-02 

1.21E-03 

6.81E-04 


Table 13: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 

2 ) 
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Here is the list of singular values: 


5.0854727e + 03 
3.8351638e + 03 
1.1897171c + 03 
5.3147287e + 02 
4.8395378e + 02 
2.6610648e + 02 
2.2399671e + 02 
1.3298740c + 02 
2.4897525e + 01 
2.1671352e + 01 
4.0442263e + 00 
1.1685319c + 00 
1.1162751c + 00 
1.0360077e + 00 
1.0219474c + 00 
1.0204436c + 00 
9.6435485e - 01 
1.8047916e - 01 
1.4189555e - 01 
5.7273237e - 02 
5.2015201e - 02 
1.8247701c - 02 
7.1876819e - 03 
4.4370659e - 03 
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Appendix 4: Complete Data for Estimation Using the 
Best Combination of Two - Parameters Initialized to Zero 
(Case 3) 

Below presents 1) a complete set of figures illustrating the estimation process, 2) a ta- 
ble describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters, and 3) a complete list of the singular values of A. 



Figure 72: Measurement Estimation - Sensed Input 1, ZWF36 (Case 3) 
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Figure 73: Measurement Estimation - Sensed Input 2, ZA8 (Case 3) 



Figure 74: Measurement Estimation - Sensed Input 3, ZA16 (Case 3) 
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Figure 75: Measurement Estimation - Efficiency and Flow Scalar 1, ZSW2 (Case 3) 





Figure 76: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 3) 
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* tO* 


Figure 77: Measurement Estimation- Efficiency and Flow Scalar 3, ZSW7D (Case 3) 
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Figure 79: Measurement Estimation- Efficiency and Flow Scalar 5, ZSW27 (Case 3) 
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Figure 80: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case 3 
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Figure 81: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 3) 



Figure 82: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 3) 
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Figure 83: Measurement Estimation - Efficiency and Flow Scalar 9, ZSW49 (Case 3) 
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Figure 85: Measurement Estimation - Bleed 1, ZWB3 (Case 3) 
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Figure 86: Measurement Estimation - Power Extraction 1, ZPWXH (Case 3) 
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Parameter 


6 


ZSW2 


SEDM2 


ZSW7D 


SEDM7D 


ZSVV41 


ZSE41 


ZSW49 


9 


ZPWXH 


Actual Value 

Mean:KF 

Mean:IF 

Variancce:KF 

Variance:IF 

3.00E-03 

6.25E-03 

6.57E-03 

1.54E-02 

1.52E-02 

0.00E4-00 

3.87E-02 

3.90E-02 

1.71E-02 

1.70E-02 

O.OOE+OO 

1.12E-01 

1.12E-01 

9.06E-03 

9.10E-03 

3.00E-03 

2.94E-03 

2.94E-03 

1.56E-05 

1.56E-05 

3.00E-03 i 

2.90E-03 

2.90E-03 

6.50E-06 

652E-06 

3.00E-03 

3.21E-03 

3.21E-03 

4.55E-05 

4.56E-05 

3.00E-03 

9.20E-04 

9.21E-04 

9.38E-05 

9.42E-05 

3.00E-03 

2.92E-03 

2.92E-03 

1.94E-06 

1.95E-06 

3.00E-03 

3.11E-03 

3.11E-03 

2.33E-06 

2.38E-06 

3.00E-03 

2.98E-03 

2.98E-03 

1.30E-06 

1.29E-06 

3.00E-03 

3.05E-03 

3.05E-03 

1.26E-06 

1.29E-06 

3.00E-03 

3.01E-03 

3.01E-03 

6.80E-06 

7.03E-06 

3.00E-03 

3.05E-03 

3.05E-03 

1.41E-06 

1.48E-06 

0.00E+00 

-7.85E-04 

-8.02E-04 

5.15E-04 

5.18E-04 

0.00E+00 

1.08E-02 

1.08E-02 

1.75E-04 

1.76E-04 


Table 14: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 

3} 
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Singular values: 


‘ 1.8399722c + 

03 ' 

1.0276994e + 

03 

5.65601 75e + 

02 

3.3739302e + 

02 

1.3176407e + 

02 

1.0460710e + 

02 

6.1757593e + 

01 

6. 554538 le + 

00 

1.4091978c + 

00 

9.3456088e - 

01 

9.2003604e - 

01 

9.0724849e - 

01 

9.0596324e - 

01 

6.0473198e - 

01 

2.1403897e - 

01 

1.4343693e - 

02 

7.0246684e - 

03 

3.9373848e - 

03 

3.0239887e - 

03 

2.6316146e - 

03 

. 1.2291035e - 

03 . 
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Appendix 5: Complete Data for Estimation Using Single 
Operating Point - Parameters Initialized to Zero (Case 4) 


Below presents 1) a complete set of figures illustrating the estimation process, 2) a ta- 
ble describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters, and 3) a complete list of the singular values of A. 



Figure 87: Measurement Estimation - Sensed Input 1, ZWF36 (Case 4) 


NASA/CR— 2006-2 14244 


207 





T»»»M ,,,* 


Figure 88: Measurement Estimation - Sensed Input 2, ZA8 (Case 4) 
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Figure 89: Measurement Estimation - Sensed Input 3, ZA16 (Case 4) 
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Figure 91: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 4) 
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Figure 93: Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D (Case 4) 
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Figure 95: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case 4) 
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Figure 96: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 4) 
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Figure 99: Measurement Estimation - Efficiency and Flow Scalar 10, ZSE49 (Case 4) 
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Figure 100: Measurement Estimation - Bleed 1, ZWB3 (Case 4) 
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Parameter 

Actual Value 

Mean:KF 

Mean: IF 

Variancce:KF 

Variance: IF 

ZWF36 

3.00E-03 

7.26E-05 

7.26E-05 

2.28E-08 

2.30E-08 

ZA8 

0.00E+00 

-2.70E-03 

-2.70E-03 

1.99E-05 

2.00E-05 

ZA16 

O.OOE+OO 

1.23E-03 

1.23E-03 

6.08E-06 

6.1 IE-06 

KawNf 

3.00E-03 



3.07E-06 

2.98E-06 

SEDM2 

3.00E-03 

2.94E-03 

2.94E-03 

4.06E-06 

4.05E-06 

ZSW7D 

3.00E-03 

2.94E-03 

2.94E-03 

6.32E-06 

6.33E-06 

llsismiiidi 

3.00E-03 



5.02E-05 

5.05E-05 


3.00E-03 

3.15E-03 

3.15E-03 

4.90E-06 

4.97E-06 

SEDM27 

3.00E-03 

3.45E-03 

3.45E-03 

4.64E-07 

4.67E-07 


3.00E-03 

3.05E-03 

3.05E-03 

5.86E-07 

5.77E-07 


3.00E-03 

3.42E-03 

3.42E-03 

1.92E-06 

1.94E-06 

ZSW49 

3.00E-03 

L83E-03 

1.83E-03 

8.08E-07 

8.78E-07 

iBMaagEM 

3.00E-03 

4.65E-04 

4.64E-04 

5.23E-07 

5.59E-07 

ZWB3 

0.00E+00 

-5.57E-04 

-5.57E-04 

6.73E-07 

6.78E-07 

ZPWXH 

0.00E+00 

-1.47E-06 

-1.47E-06 

6.08E-12 

6.14E-12 


Table 15: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 
4) 


Singular values: 

' 3.0135163e + 03 " 
1.9332938c + 03 
6.8244808e + 02 
1.7598337e + 02 
1.0068188c + 02 
9.4230100e + 00 
7.5684915e + 00 
3.5334475e + 00 
1. 03204 13e + 00 
9.4464178e - 01 
8.0679049e - 01 


(6.4) 
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Appendix 6: Complete Data for Estimation Using the 
Best Combination of Three - Parameters Initialized to 
0.5 (Case 5) 

Below presents 1) a complete set of figures illustrating the estimation process, and 2) a 
table describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters. 
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Figure 102: Measurement Estimation - Sensed Input 1, ZWF36 (Case 5) 
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Figure 103: Measurement Estimation - Sensed Input 2, ZA8 (Case 5) 



Figure 104: Measurement Estimation - Sensed Input 3, ZA16 (Case 5) 


NASA/CR— 2006-2 14244 


218 






Figure 105: Measurement Estimation - Efficiency and Flow Scalar 1, ZSW2 (Case 5) 



Figure 106: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 5) 
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Figure 109: Measurement Estimation- Efficiency and Flow Scalar 5, ZSW27 (Case 5) 
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Figure 112: Measurement Estimation - Efficiency and Flow Scalar 8, ZSE41 (Case 5) 
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Figure 113: Measurement Estimation - Efficiency and Flow Scalar 9, ZSW49 (Case 5) 



Figure 114: Measurement Estimation - Efficiency and Flow Scalar 10, ZSE49 (Case 5) 
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Figure 115: Measurement Estimation - Bleed 1, ZWB3 (Case 5) 
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Figure 116: Measurement Estimation - Power Extraction 1, ZPWXH (Case 5) 
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SEDM7D 


ZSW27 


SEDM27 


ZSW41 


ZSE41 


ZSW49 


ZSE49 


ZWB3 


Actual Value 


3.00E-03 


O.OOE+OO 


O.OOE+OO 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


3.00E-03 


O.OOE+OO 


O.OOE+OO 


MeanrKF 

Mean:IF 


-1.98E-02 

-1.94E-02T 

iBBHasamaiai 

2.95E-03 

2.95E-03 

2.99E-03 

2.99E-03 

2.97E-03 

2.96E-03 

2.84E-03 

2.81E-03 

3.00E-03 

3.01E-03 

2.96E-03 

2.97E-03 

3.00E-03 

3.00E-03 

3.10E-03 

3.10E-03 

2.99E-03 

3.00E-03 

3.00E-03 

3.01E-03 

1.64E-04 

1.49E-04 



VariancerKF 


1.33E-02 


1.08E-03 


2.69E-03 


2.36E-07 


1.13E-07 


7.38E-07 


9.75E-07 


4.37E-08 


7.89E-08 


8.18E-08 


1.56E-07 


1.80E-07 


1.70E-07 


2.36E-05 


2.85E-02 


Table 16: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 

■5) 
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Appendix 7: Complete Data for Estimation Using the 
Best Combination of Two - Parameters Initialized to 0.5 
(Case 6) 

Below presents 1) a complete set of figures illustrating the estimation process, and 2) a 
table describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters. 
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Figure 117: Measurement Estimation - Sensed Input 1, ZWF36 (Case 5) 
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Figure 118: Measurement Estimation - Sensed Input 2, ZA8 (Case 6) 
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Figure 119: Measurement Estimation - Sensed Input 3, ZA16 (Case 6) 
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Figure 120: Measurement Estimation - Efficiency and Flow Scalar 1, ZSW2 (Case 6) 



Figure 121: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 6) 
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Figure 122: Measurement Estimation- Efficiency and Flow Scalar 3, ZSW7D (Case 6) 



Figure 123: Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D (Case 6) 


NASA/CR— 2006-2 14244 


230 






Figure 124: Measurement Estimation- Efficiency and Flow Scalar 5, ZSW27 (Case 6) 
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Figure 125: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case ( 
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Figure 126: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 6) 
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Figure 128: Measurement Estimation - Efficiency and Flow Scalar 9, ZSW49 (Case 6) 



Figure 129: Measurement Estimation - Efficiency and Flow Scalar 10, ZSE49 (Case 6) 


NASA/CR— 2006-2 14244 


233 





0 OS 


1 S 2 2.5 3 

(Me) 


iJ 


* JO* 


Figure 130: Measurement Estimation - Bleed 1, ZWB3 (Case 6) 
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Figure 131: Measurement Estimation - Power Extraction 1, ZPWXH (Case 6) 
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Parameter 


ZWF36 


ZA8 


ZA16 


ZSW2 


SEDM2 


ZSW7D 


SEDM7D 


ZSW27 


SEDM27 


ZSW41 


ZSE41 


ZSW49 


ZSE49 


Z\VB3 


ZPWXH 


Actual Value 

MeanrKF 

Mean:IF 

VarianccerKF 

Variance:IF 

3.00E-03 

1.15E-01 

1.15E-01 

1.67E-02 

1.65E-02 

O.OOE+OO 

1.53E-01 

1.53E-01 

2.82E-02 

2.83E-02 

0.00E+00 

2.74E-01 

2.74E-01 

2.24E-02 

2.25E-02 

3.00E-03 

3.93E-03 

3.93E-03 

3.03E-05 

3.04E-05 

3.00E-03 

3.30E-03 

KUTiTalfil 

9.47E-06 

9.50E-06 

3.00E-03 

1.78E-03 

1.77E-03 

5.12E-05 

5.16E-05 

3.00E-03 

5.41E-03 

5.45E-03 

9.66E-05 

9.73E-05 

3.00E-03 

3.39E-03 

3.39E-03 


6.05E-06 

3.00E-03 

' 2.92E-03 

2.91E-03 

2.55E-06 

2.60E-06 

3.00E-03 

wsBssmmmi 

1.52E-06 

1.51E-06 

3.G0E-03 

3.06E-03 

3.06E-03 

3.59E-06 

3.64E-06 

3.00E-03 


1.24E-05 

1.26E-05 

3.00E-03 

3.12E-03 

3.12E-03 

1.87E-06 

1.93E-06 

0.00E+00 

3.23E-03 

3.17E-03 

1.64E-03 

1.65E-03 

0.00E+00 

4.67E-01 

4.67E-01 

1.74E-03 

1.75E-03 


T;i!ile 17: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 

f> * 
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Appendix 8: Complete Data for Estimation Using Single 
Operating Point - Parameters Initialized to 0.5 (Case 7) 


Below presents 1) a complete set of figures illustrating the estimation process, and 2) a 
table describing the mean and variances of the parameter estimates of both the Kalman and 
Information filters. 


mV Waw m i m to* t»i 1 uMnf on* oMranng pen 



Figure 132: Measurement Estimation - Sensed Input 1, ZWF36 (Case 7) 
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Figure 133: Measurement Estimation - Sensed Input 2, ZA8 (Case 7) 
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Figure 134: Measurement Estimation - Sensed Input 3, ZA16 (Case 7) 
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Figure 135: Measurement Estimation - Efficiency and Flow Scalar 1, ZSW2 (Case 7) 



Figure 136: Measurement Estimation - Efficiency and Flow Scalar 2, SEDM2 (Case 7) 
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Figure 138: Measurement Estimation - Efficiency and Flow Scalar 4, SEDM7D (Case 7) 
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Figure 140: Measurement Estimation - Efficiency and Flow Scalar 6, SEDM27 (Case 7) 
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Figure 141: Measurement Estimation - Efficiency and Flow Scalar 7, ZSW41 (Case 7) 



Figure 142: Measurement Estimation - Efficiency and Flow Scalar 8, ZSE41 (Case 7) 
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Figure 144: Measurement Estimation - Efficiency and Flow Scalar 10, ZSE49 (Case 7) 


NASA/CR— 2006-2 14244 


243 






Kaanan ana mtormanon P*ar» tor mpu 14 uamg on* o pr a itag po« 



*10* 


Figure 145: Measurement Estimation - Bleed 1, ZWB3 (Case 7) 


Kaanan ana ta «ow wH w» p»«ra tar »tni 1 S uaaig on* op*r a nng port 



Figure 146: Measurement Estimation - Power Extraction 1, ZPWXH (Case 7) 
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Parameter 


ZWF36 


ZA16 


ZSW2 


SEDM2 


ZSW7D 


SEDM7D 


ZSW27 


SEDM27 


ZSW41 


ZSE41 


ZSW49 


ZSE49 


'B3 


ZPWXH 


Actual Value 

Mean:KF 

Mean:IF 

Variance:KF 

Variance: IF 

3.00E-03 

4.94E-01 

4.94E-01 

2.11E-08 

wmmmm 

O.OOE+OO 

7.66E-01 

7.66E-01 

1.85E-05 

1.86E-05 

0.00E+00 

3.98E-01 

3.98E-01 

5.68E-06 

5.70E-06 

3.00E-03 

2.09E-02 

2.09E-02 

3.34E-06 

3.26E-06 

3.00E-03 

4.03E-03 

4.03E-03 

4.61E-06 

4.61E-06 

3.00E-03 

7.02E-03 

7.02E-03 

6.81E-06 

6.83E-06 

3.00E-03 

1.14E-01 

1.14E-01 

4.65E-05 

4.68E-05 

3.00E-03 

4.42E-02 

4.42E-02 

5.40E-06 

5.46E-06 

3.00E-03 

-1.55E-01 

-1.55E-01 

5.23EMJ7 

5.27E-07 

3.00E-03 

1.32E-02 

1.32E-02 

5.78E-07 

5.64E-07 

3.00E-03 

-3.99E-02 

-3.99E-02 

1.89E-06 

1.91E-06 

3.00E-03 

2.22E^01 

2.22E-01 

8.28E-07 

8.95E-07 

3.00E-03 

4.22E-01 

4.22E-01 

5.32E-07 

5.66E-07 

0.00E+00 

5.61E-01 

5.61E-01 

6.24E-07 

6.27E-07 

0.00E+00 



5.63E-12 

5.67E-12 


Table 18: Mean and Variance of Parameter Estimates of Kalman and Information Filters (Case 

7) 
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Introduction 


The main purpose of the work reported here is to design a controller or tune a PE) 
controller for the GE90 CPFN fuel system. The secondary purpose is to analyze the 
stability and performance of the controlled system. 

The model for the fuel system that was used for this work was supplied by the Controls 
Center of Excellence Department of General Electric’s Aircraft Engines Engineering 
Division. It is an Excel model and it contains some basic fluid compression and fluid 
momentum effects. 

The fuel system model is partially described in the next section followed by a section 
listing the controller goals. The section after that discusses the basic controller that was 
included in the spreadsheet as well as a new controller that is proposed in this paper. The 
stability of the fuel system is then detailed, followed by simulation results showing the 
performance of the controlled system with the basic and proposed controllers. 

Fuel System Model 

Components 

The components of the fuel system are the main pump, boost pump, bypass valve, fuel 
metering valve (FMV), pressurizing valve (PV), controlled pressure fuel nozzle (CPFN) 
pressurizing servovalve (PSV), and CPFN valves. The purposes of these are described 
below. 

The main pump should be capable of delivering enough fuel for the combustor at 
maximum usage and also deliver fuel to the PV and CPFN PSV for control purposes. The 
goal of the bypass valve is to bypass pump flow back toward the pump inlet to maintain a 
constant prescribed pressure drop across the FMV. However, during transients, the 
pressure drop may deviate by about 5 psi. The flow through the FMV is simply 
proportional to the flow number of the FMV when the pressure drop is at the prescribed 
value. The purpose of the CPFN PSV is to supply the signal side of the CPFN valve with 
the proper pressure and flow to move the CPFN piston to follow a demanded position. 
The purpose of the PV is to keep the boost pressure at least 250 psi above the supply 
pressure. 
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Equations in EXCEL 

The following equations were extracted from the EXCEL model. They were included 
here mostly to document the model and can be skipped without loosing the main focus of 
this paper. 

Engine: 

P 36 = ((607-90.9) / (24639-2440))*WF 36 + 54 

XN 2 = ((9 1 97 -6066)/(24639-2440)* WF 3 6 + 4700)*6498 1/5722 

PSV: 


Pc = Pboost + 160 

DPpv = IF(MODEpsv > l,Pc - PiguessiO) 

V PV = IF(DPp V > 0,77*K 2 7 *SQRT(DPpv)/K25,-77*K27*SQRT(-DPpv)/K 2 5)) 

DVDT = (Vp V - Vpv(k-l)) / T 

F NEr = DVDT*K 28 / (32.12*12) 

MODE = (1 for stuck, 2 for forward, 3 for backward) 

FORCE = (P 2 - P c + DPp V )*K25 - K 22 - K 2 i*X PV (k-l) 

FRICTION = IF MODE=l THEN max(-K 29 ,(min(K 29 ,FORCE))) 

IF MODE=2 THEN -K 30 
IF MODE=3 THEN +K 30 

ERRXpv = FORCE + FRICTION - F NE t 

V pv = IF(MODEpsv < 2,0,V PV ) 

Xpv = max(K 2 4 ,min(K 23 ,Xpv(k- 1 ) + Vpv 2 -.02))) 2 + P 22 

FNpv = (16689*0.6*(2.14*X PV 2 - 0.02)) 

P 2 = P 22 + (W fm /FNpv ) 2 
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CPFN: 


Fsprdmg = K*Xcpfn + PRELOAD 

FRICTION = IF MODE=l THEN max(-FST,(min(FST,FORCE))) 

IF MODE=2 THEN - FSL 
IF MODE=3 THEN +FSL 

Fnet = FORCE + FRICTION + Fsrping 

A — Fnet / MASS 

V = V(k-1) +a*T 

Xcpfn = XcPFN(k-l) + V*T 

XdMD.CPFN = f(XoMD,FMv) 

A5S4 = F 19 *X SENS ,cPFN(k-l) + H 19 *XcPFN(k-l) 

XsENS = A5S4 + Il 9 *XcPFN 

ERRXcpfn = Xsens.cpfn - Xdmd.cpfn 

ERRX2cpfn = ERRXcpFN(k-l) - max(min(ERRXcpFN(k-l), 0.005,-0.005) 
ERRX3cpfn = max(min(ERRX2cpFN>0.01,-0.01)) 

STBIcpfn = STBlcpFN(k-l) + KIcpfn*ERRX3cpfn 


STB2cpfn = STB1 - Hu 
INTPART = STB2cpfn 

KPpsv = IF(TMC4 < 0,E 2 5*sqrt(100/(Pi(k-l)-20)),E25*sqrt(I00/(P s (k-l)-P 1 (k-l))) 

PROPPART = KPpsv *ERRXcpfn 

TMC4 = min(100,max(-100,PROPPART + INTPART)) 

TMC4N = TMC4 + H„ 

TMC4H = TMC4H(k-l) + min(TMC4N - TMC4H(k-l) + 0.5*HYS C pfn,0) 

+ max(TMC4N - TMC4H(k-l) - 0.5*HYS C pfn,0) 
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Pi = P 2 + 0.0278*K14* dQj/dt / ((386*7t*K15 A 2)/4 + sign(Q P sv)*0.001892*Q P sv A 2) 

Qexpi = H 13 *H 14 *SQRT(abs((P 1 -20)/100))*TMC4H + (l-H 13 )*Q P sv(K-l) 

Qexp2 = Hi 3 *H 14 *SQRT(abs((P s -Pi)/100))*TMC4H + (1-H 13 )*Q PSV (K-1) 

Qexp3 = -Hi 3 *H 14 *SQRT(abs((Pi-20)/100))*TMC4H + (1-H 13 )*Q PSV (K-1) 

Qexp4 = -Hi 3 *H 14 *SQRT (abs((Ps-P 1 )/ 1 00))*TMC4H + (1-H 13 )*Q PSV (K-1) 

IF (Pi > 20) 

THEN 

Qbias = -0.0465*SQRT(P r 20) 

IF Qexpi < 0 
THEN Q P sv = Qexpi 
ELSE Q P sv = Qexp2 

ELSE 

Qbias = 0.465*SQRT(20-P,) 

IF Qexp3 < 0 
THEN Qi>sv = Qexp3 
ELSE Q P sv = Qexp4 

END 

Qpsv = Qpsv + Qbias 
q2v = from solver 
Q 2 = 30* q2v 

dP 2 / dt = (Qpsv - Q 2 )*75000/(K 14 *7 i*(Ki 5 a 2)/4) 

P 2 = P 2 (k-l)+T* dP 2 /dt 

dQi / dt = (Qpsv - Qpsv(k-1))/T 

Physical Equations for the CPFN Valve 

The force balance equation for the piston position is 

(P T - Pb)A = F f + (Kx + preload) + ma . 

where: 

P T = pressure at the top of the piston (manifold side) 

P B = pressure at the bottom of the piston (signal side) 
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A = area 

Fp = force of friction 
K = spring constant 
x = piston position 
preload = spring force at x=0 
m = mass 
a = acceleration 

In the limit, as the mass of the piston goes to zero, the term “ma” goes to zero. The 
equation then becomes 

PdiffA = Fp + (Kx + preload) . 

where: 


Pdiff — Pt -Pb 


The flow through the hole in the piston is 
Qt = FN*sqrt(P D(F F) 

Moving Piston 

When the piston is moving, the velocity of the piston is related to the net flow in the 
bottom half of the cylinder. Since the bulk modulous of liquids, like jet fuel, is high, a 
good approximation can be obtained by neglecting fluid compression: 

AV = Q b -Qt or 

AV = Q b -FN*sqrt((Kx + preload + F f )/A) . 

where: 

V = velocity 

Qt = flow from the manifold side (TOP) of the CPFN valve to the signal side through the 
piston orifice 

Qb = flow out of the signal side (BOTTOM) of the CPFN piston towards the CPFN PSV 
When the piston has positive velocity, F f = +Fst, so the equation relating Qb to velocity is 

V = Qb/A - FN/A* sqrt((Kx + preload + F S t)/A) . 

When the velocity is negative, 

V = Qb/A - FN/A* sqrt((Kx + preload - F S t)/A) . 
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Therefore, velocity is proportional to the flow out of the bottom of the piston with a bias 
term that is a function of the position. When x is large enough, the orifice flow is not 
affected much by the direction of the velocity. 

Plots of the above simplified equations were made in EXCEL and it was shown that they 
are very good approximations to the EXCEL equations with fluid compression and 
nonzero mass. 

Static Piston 


When the piston is static, 

Qt = FN*sqrt((Fp + Kx + preload)/A) 

where the frictional force in the force balance equation satisfies the inequality 
-Fst < F f <Fst • 

From the equation for Qt, it is seen that for a particular position, the range of Ff implies a 
range of allowable flows. 

CPFN Transfer Functions 


A linear model of the CPFN loop was calculated to determine the linear stability. It was 
calculated using the physical equations for the CPFN piston (of the previous section), the 
Excel equations for the sensor and actuator, and the control laws that apply when the 
piston is moving. 

To determine stability of a linear model, the external inputs to the loop are not relevant. 
The stability is just a function of the open loop gain. The transfer functions, by definition 
are the ratio of the transform of the output of a block divided by the transform of the 
input, neglecting effects of the output on external inputs. 

The transfer function from CPFN velocity to position is 

Gi(z) = X(z)/V(z) = T/(z-l). 

where: 

T = integration time step 

The transfer function of the sensor is 

G 2 (z) = Xsens(z) / X(z) = (1 - 0.6065) / (z - 0.6065) . 
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It therefore has a steady-state gain of unity. The transfer function of the controller when 
the piston is moving is 

G 3 (z) = TMC4(z) / X SE ns(z)= ( K p *K p ,psv z + (K,*T-K P *K P>P sv) ) / (z-1) . 
where: 

K P = proportional gain when piston is moving 
Ki = integral gain when piston is moving 
K PjP sv = gain that is a function of Pi and Ps 
TMC4 = CPFN torque motor current 

The input to the controller is the sensed position and the output is the torque motor 
current. The gain from the CPFN torque motor current, TMC4, to the flow, Q P sv is 

G 4 (z) = Qfsv(z) / TMC4(z) = g PSV *0.18*(l - 0.6065) / (z-0.6065) . 

where: 

gpsv = gain for the torque motor current, this is a function of P] and P s 

The steady state gain of G 4 (z) is gpsv*0.18. From the previous section, the equation 
relating velocity to flow (neglecting fluid compression) is AV = Qb -Qt- Since there are 
30 CPFN valves 


Qb= -Qpsv/30 


and the equation becomes 


AV = - Qrev /30 - Qt- 


For the linear model, Qt is treated as a disturbance. Therefore, the transfer function from 
the flow to velocity is 

G5(z) = V(z) / Qpsv(z) = -1 / (30* A). 

The factor KPpsv in the controller was designed to nearly cancel out the effect of the 
variable gain, g PS v, from the TMC to the signal flow. If the TMC were not limited, and if 
the factor KP PS v multiplied both the proportional and integral terms, KPi>sv would 
exactly cancel out the gain change from the TMC to the signal flow. 
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Controller Goals And Issues 


Controller Goals 

1. Control the CPFN piston to make single stage change within 0.5 seconds. 

2. The steady-state error of the CPFN piston should be small. 

3. When stage changes are not requested, the controller should keep the CPFN piston 
still, if possible. If it does not hold the piston still, the movements should either be 
small or driven back to the set point quickly. This applies particularly to movements 
due to changes in the fuel manifold pressure. 

4. Keep dips in the combustor fuel flow due to CPFN piston flow small. 

5. The controller should have good stability margins with the linear model of the fuel 
system. 

6. The fuel nozzle pressure drop should be greater than 20 PSI. 

Uncommanded Piston Movements 


When the commanded value of the CPFN piston is held constant, the piston, 
unfortunately can still move due to changes in the manifold pressure. This can occur 
when the demanded fuel flow is a ramp and the demanded CPFN piston position is a 
series of steps. The position sensor will not indicate the buildup of pressure because the 
position will not change until the static friction is exceeded. Therefore, advance notice of 
the movement, if desired, will have to come from pressures. 

Controllers 

Basic Controller 

Analex received the spreadsheet with a proportional and integral controller with position 
error as input. The position error is the demanded position minus the sensed position. For 
the integral term, this error was put through a dead band of 0.005 inches, multiplied by 
the integral gain (Ki), integrated, and then the output was rate limited to a change of 
±0.01 per integration time step (T) of 0.01 seconds. 

For the proportional part, the error was multiplied by a gain that is a function of Pi and P s . 
This gain, KP P sv, was chosen to be the reciprocal of the gain from the torque motor 
current (TMC4) to the pressurizing servovalve current, Q PSV . 
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Proposed Controller 

This report proposes to use the previously mention controller when the piston is moving 
and to use a controller based on the estimated value of the CPFN static friction when the 
piston is static. (The actual value of the friction is known in the simulation, but with the 
physical system, it must be estimated). The reason for this, is that a controller sensing 
only position, can not determine when the pressure difference across the cylinder is 
building up and is about to cause unwanted piston movement. By monitoring the 
estimated friction, it is possible to predict when the piston is going to break away from its 
position. If the direction of movement of the friction is not known, then it makes sense to 
attempt to keep the friction near zero because this will allow the largest deviation without 
breakaway. If the direction in which the friction is likely to move is known, then the 
target friction should be set to be near the static friction limit in the opposite direction. 

According to a previous section, when the piston is static, the frictional force should be 

F st = A*(Pman - PI) - (K*X - pre-load). 

Therefore, driving the following error signal to zero should drive the static friction to zero 
also. 


errorop = Pman,sens - PI sens - (Kcont*Xsens - pre-loadcoNT) / Acont • 


where: 

errorop = error signal for the pressure feedback controller (delta pressure controller) 
Pman,sens = sensed value of P M an 
PI sens = sensed value of PI 

Kcont = value of the spring constant used by the controller 
pre-loadcoNT = value of the pre-load used by the controller 
Acont = value of the area used by the controller 

Since the position feedback controller works moderately well alone, i.e. for a dynamic 
and static piston, it is proposed that the pressure feedback controller be used only when 
the measurements provide a very high probability that the piston is in the static mode. 
The following criteria are recommended to be used to determine when to activate the 
pressure feedback controller (static controller): 

1 position error below a given threshold 

2 velocity has been below a threshold for “n” consecutive time periods or 
low passed filtered velocity signal is below a threshold 

It is recommended that both of these be satisfied when using pressure feedback. The 
velocity criteria may introduce a small delay in the turning on the pressure feedback 
controller, but this will not cause problems. There will be no delay switching from the 


NASA/CR— 2006-2 14244 


258 



static controller to the dynamic controller because even though criteria #2 may still be 
satisfied, criteria #1 will not be and the position feedback controller will be turned on. 

To have the output of the integrator properly initialized when the pressure feedback 
controller is turned on, the same integrator can be used for both controllers. However, it is 
fine to use different dead bands and rate limiters. 

Since the sensed pressures and the spring parameters will be slightly different than the 
actual values, the actual friction force will not be exactly zero when the feedback signal is 
zero. However, a full dispersion analysis is beyond the scope of the present task. 

Flow Control 


Another method to drive the frictional force to zero is to control the flow out of the signal 

end of the cylinders. This could be accomplished as follows: 

1 calculate a pressure drop across the piston, based on the position, to cancel out the 
spring force (i.e. no static friction) 

2 calculate the flow through the orifice corresponding to this pressure drop, denote this 
as the demanded flow (this depends on an accurate value for the orifice flow number) 

3 calculate the actual flow by measuring a pressure drop in the path 

4 drive a simple controller with the error between demanded flow, #2, and the actual 
flow, #3 

The method in the previous section was chosen over this method because this method 

relies on an accurate value for the orifice flow number and on measuring the flow. 

Fuel Demand Modification 


When the CPFN piston changes stages a flow drop to the combustion chamber occurs 
which is equal to the area times the velocity. The FMV does not have the bandwidth to 
cancel out this disturbance. In fact, even if the demand is increased as an attempt cancel 
this disturbance, the disturbance is still not canceled, but it can be significantly reduced. 
The method that was employed for this work was to modify the fuel demand by adding 
half a period of a sine wave to the demanded fuel flow. To change the demand for a time 
interval of W seconds, the following equation was used. 

WF DM d,add = A*sin(co*(t - (t stage - W/2))) for (t stage - W/2) < t < (t stage + W/2) 
where: 

tstage = anticipated time for the CPFN piston to start a stage change 

T=2*W 

f=l/T 

co = 2*n*f 
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Reasonable values for the parameters are A=800 lbm/hr and W=0.3 seconds. A square 
pulse could also be used, but it would change the demand abruptly. 

The main difficulty with this modification is that the stage change must either be known 
in advance or be delayed by a time period equal to W. And since the demanded stage may 
be a function of the actual fuel flow, changing the demanded flow will change the actual 
flow, which may change the demanded and actual time of the stage change. 

Stability of the Linear Model 

Nonlinear stability is often difficult or even intractable for problems with nonlinear 
blocks that have equations that change abruptly in time, as does the CPFN pressure 
regulator. Because of this, it was decided to determine the stability of linearized versions 
of the nonlinear system (at several operating points) and also to investigate the stability of 
the nonlinear system in the time domain with simulations. 

To determine the stability, only the CPFN loop was included in the model. This loop 
consists of the CPFN regulator and the CPFN valve. The nonlinearities of this loop are in 
the equations for Qpsv and KPpsv- The equations for these variables are functions of 
sqrt(100/(Pi(k-l)-20)) or sqrt(100/(Ps(k-l)-Pi(k-l))) depending on the direction of the 
main part of the flow of Qpsv- The main part does not include the bias flow. For each of 
the three parameters: flow direction, Pi, and Ps there is a linear model. This linear model 
is the set of transfer functions given in the section “CPFN Transfer Functions”. The 
manifold pressure is treated as a disturbance in these models and does not affect the 
stability. 

The model was reparametized to be a function of flow direction, Pj, and the difference 
between the pressures “dp” where dp=Ps-Pi. Ranges for Pi and dp were chosen and the 
stability was determined for the eight cases that are obtained when each parameter can be 
at one of its extremes. The extremes that were chosen for the stability analysis are 40 and 
200 for Pi, 300 and 500 for dp, and forward or backward for the main part of the flow 
Qpsv- 

The linear models were analyzed in the frequency domain and with root locus. The roots 
were stable for all eight cases and the stability margins for a controller with Kp=2200 and 
Ki=24 were in the ranges listed below. 

2.43 < Gain Margin < 2.44 

7.74 < Gain Margin (db) < 7.75 

41.7 < Phase Margin < 41.8 deg 
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The ranges are very tight, because the controller gain at these frequencies is mostly based 
on the proportional term, and the factor KPpsv in this term nearly cancels out the effect of 
the changing gain from the TMC to the signal flow. 

Typical gain and phase margin goals are 6 db and 30 deg. The calculated gain margin for 
this fuel system implies that the system would be stable for plant a uncertainty factor of 
2.43. The phase and gain margins occur at frequencies of 15.4 and 30.3 rad/sec 
respectively. Therefore, this phase margin shows that the linear system will be stable at 
2.45 hz with an unmodeled delay of 0.047 seconds. 

The stability of the linear models was determined by the Matlab function in Appendix A. 

Results 

Position Control 

Results of a simulation using the basic proportional plus integral position controller with 
gains KP=2200 and KI=24 are shown in Appendix B. Figure B-la shows the demanded 
and actual CPFN piston position as a function of time. It changes from stage 1 to 2 to 3. 
Figure B-lb shows that the velocity becomes nonzero around 0.25 seconds. This is due to 
an increase in the pressure difference across the cylinder, which is due to an increase in 
manifold pressure that is caused by an increase in demanded fuel flow. The demanded 
fuel flow is a ramp function in this simulation. Similarly, the piston breaks loose around 
1.05 and 1.7 seconds. Since the feedback signal to the controller is position, the best it 
can do is to react after the piston has moved. The next section will illustrate the 
performance of a controller, to be used when the piston is static, that senses a pressure 
difference and drives it towards zero to reduce the number of breakaways. 

The model was received by Analex with the CPFN piston mass set equal to 1.0 lbm. 
Current specifications call for the mass to be 0.138 lbm. When the piston is still, the mass 
is irrelevant. When the piston is in motion, the pressure, spring, and frictional forces are 
much larger than the inertial force of the piston (because of its small mass) and a very 
good approximation to the transient response can be obtained by neglecting the inertial 
force. This was verified by performing three simulations in which the CPFN piston mass 
was set to the values of 0.138, 1.0, and 10 lbm. The plots were indistinguishable to the 
human eye. Since the EXCEL solver frequently had difficulty with a piston mass of 
0.138, a value of 1.0 was used in the runs described in this paper. 

Position and Pressure Control 

Results of a simulation using the position feedback controller while the piston is moving 
and a pressure feedback controller while the piston is static are shown in Appendix C. 
Figure C-la shows that the piston moves from stages 1 to 6. Figure C-lb depicts the 
velocity. After the first stage transition the velocity does not reach zero so the static 
controller was not activated. It was activated after the 2 nd , 3 rd , and 4 th transitions and it 
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prevented the piston from breaking loose. The breakaway that occurs after the 5 th 
transition, happens before the static controller is turned on. For this simulation the static 
controller was not allowed to become active until a lowpassed velocity signal was within 
certain limits. The breakaway that occurs before the 1 st transition is due to a change in the 
dynamics of the PSV flow. In the model, this is manifested as a change in the PSV flow 
equation. The static controller was not predicted to prevent a breakaway due to this 
change in dynamics which occurs infrequently. The actual value of the friction in the 
simulation (as opposed to the estimated value based on pressure transducers) is shown in 
Figure C-lc. It is seen that the static controller drives the friction to zero when is it active 
(zero velocity). It is also seen that the friction oscillates much more at higher stages. The 
oscillation can be greatly reduced by either reducing the gain at higher stages (gain 
scheduling) or by using pressure transducers with higher bandwidths. The pressure 
transducer had a cut-off frequency of 8 hz in this simulation. 

Figure C-2a shows the value of the controller switch and the position error. When the 
switch is zero the dynamic controller is active, when it is one the static controller is 
active. Figure C-2b shows the proportional and integral outputs of the active controller, 
i.e., whichever controller is active at the time. The integral part is on top. The 
proportional part is not strictly a constant times the position when the dynamic controller 
is on due to limiting and hysteresis. 

Figure C-3a shows the controller switch and the error signal used by the pressure 
controller. This error signal has a term due to the piston position and preload. This plot 
clearly shows how the error signal is driven to zero when the static controller is active. It 
also shows that the error signal is at a constant of FSL/A when the piston position is at 
stage two but moving slowly. This occurs between 0.8 and 1.0 seconds. This illustrates 
that the pressure controller error signal (and hence the static friction) does not move to or 
stay near zero without pressure feedback control. Figure C-4a depicts the FMV demanded 
fuel flow, the FMV flow, and the combustor flow (WF36). The demanded flow is 
constant until 0.1 seconds an then a ramp after that. The actual fuel flow has slight 
“bumps” (deviations from a straight line) due to the fact that the bypass filter does not 
keep the pressure drop across the FMV exactly at the target value during transients. The 
FMV flow is the highest trace at four seconds. The combustor flow has larger bumps than 
the FMV flow due to CPFN piston movements. Figure C-4b shows the error between the 
FMV demand and the combustor flow. By observing the portions between the transients, 
it is seen that the enror is decreasing with time in a smooth and predictable fashion. 
Because of this, it seems appropriate to add a correction to the demand so that the actual 
flow is much closer to the demand. The piston movement for the transition from stage 2 
to 3 causes a fuel flow error of about 900 lbm/hr. The next three transitions cause errors 
of about 800, 1000, and 1200 lbm/hr respectively. Modification of the demand to reduce 
these drops is covered in a later section. The last plot on this page, C-4c, shows the fuel 
nozzle pressure drop. It stays above the requirement of 20 psi. 

Figure C-5a is a plot of the position error. It is seen again that the error is changing 
between the first and second transitions. The error for the next four stages is within 


NASA/CR— 2006-2 14244 


262 



±0.005 inches. Plot C-5b was made to determine the usefulness of the error signal that is 
the input to the pressure feedback controller. The error signal uses a measured value of 
Pi, which is the pressure of the output of the CPFN PSV, instead of the actual value of P3, 
which is the signal pressure of the CPFN valve. If the error signal used was the actual 
value of P3, then the trace in the figure would be zero when the piston is static. However, 
it is seen that for stage 6 (4.0 to 4.5 seconds) especially, the error oscillates about zero 
with a peak magnitude of about 20 psi and a frequency of about 10 hz. What can not be 
determined from this plot alone is how much of the deviation from zero is caused by 
pressure differences from Pi to P3 and how much of the deviation is caused by sensor or 
transducer errors. The transducer is modeled as a low pass filter. The cut-off frequency 
for this simulation is 8 hz. The next section contains results from a simulation in which 
the cut-off frequency is 48 hz. 

The gains that were used for the static controller are Kp=0.5 and Ki=0.1. 

Position and Pressure Control with a High Bandwidth Pressure Transducer 

It was stated earlier that two methods could be employed to reduce the oscillations of the 
pressure feedback controller for stage six. The simulation described in this section uses a 
pressure transducer with a larger cut-off frequency, namely 48 hz. When the pressure 
feedback controller is active for stage six the oscillations of the static friction are seen in 
Figure D-lc to be much smaller than they were with the 8 hz transducer (compare to 
Figure C-lc). The controller switch and the error signal for this run are plotted in Figure 
D-3a. 


As expected, the plots of Figure D-4 are quite similar to those of C-4. By comparing the 
controller switch of Figure D-3a with the plot of Figure D-5b, it is seen that when the 
static control is active, its feedback signal is quite good. Therefore, the feedback signal is 
good for the intervals: 

0.2 to 0.4 
1.45 to 1.6 

2.1 to 2.75 

3.2 to 3.65 
3.95 to 4.7 

This demonstrates that increasing the bandwidth of the pressure transducer is a very 
effective method to reduce the oscillation at high stages. 
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Position Control with Decreasing Flow 

This section describes a simulation in which the stages change from 6 to 1 using only the 
position feedback controller. The demanded and actual CPFN piston positions are shown 
in Figure E-la. It is seen that the actual position overshoots the demanded position during 
the first and second transition. This did not happen going in the opposite direction, i.e. 
from stage 1 to stage 6. Figure E-lc shows that without the pressure feedback controller, 
the friction hits the static limit and the piston moves as shown in Figure E- lb at 1, 2.9, 
and 4.3 seconds. Figure E-2a shows that the controller was intentionally held at zero to 
illustrate the effectiveness of position only control going from stages 6 to 1. This figure 
also shows that the error is large for stage 4 (3.2 to 3.7 seconds). Figure E-2b shows that 
for this time period, the integral output is going down linearly and presumably would 
correct the position error if the demand for stage 5 did not occur. 

Figure E-4a shows the demanded and actual FMV flows as well as the combustor flow. 
Figure E-4b shows the error. The fuel nozzle pressure drop is seen in Figure E-4c to stay 
above the requirement of 20 psi. The position error is shown in Figure E-5a. 

Position Control with Fuel Demand Modification 


Results of a simulation in which the fuel demand was modified are shown in Appendix F. 
The demanded fuel flow was modified for the first stage change which is seen in Figure 
F-la to occur at about 0.4 seconds, but it was not modified for the second stage change 
that occured at about 1.2 seconds. The fuel demand modification shown in Figure F-5c is 
seen to be half a period of a sine wave as described in a previous section. The peak error 
for the transition is about 1400 lbm/hr for the stage change without the modification and 
about 900 lbm/hr with the modification (see Figure F-4b). 

Stage Change Position Errors 

The distance between ports in the CPFN cylinder is 0.103 inches, so that the total error of 
the piston with all dispersions and control system transient errors needs to be less than 
half of this value or ±0.051 inches. GE has conducted a preliminary study on the piston 
error due to dispersions such as variances in the spring constant, spring preload, orifice 
location, friction, and LVDT accuracy. The simulations described here only show the 
position errors due to the dynamics of the system with the controller and the nonlinear 
actuator. 

Table 1 lists the position errors for the simulations of Appendices C and E. In Appendix 
C, the CPFN piston was successively stepped from stage 1 to stage 6. The position and 
pressure feedback controllers were both used in this appendix. For Appendix E, only the 
position feedback controller was used, but the position errors listed in the table were 
taken before the breakaways (see Figure E-la and E-5a). In this appendix the CPFN 
piston was successively stepped from stage 6 to 1. 


NASA/CR— 2006-2 14244 


264 



Table 1 


Stage Change Position Errors 



Appendix C 

Appendix E | 

| Stage 

Time 

Position Error 
(0.001 inch) 

Time 

Position Error 
(0.001 inch) 

1 

N/A 

N/A 

5 

-6.5 

! 2 

N/A 

N/A 

4 

1.3 

3 

1.5 

1.9 

3.5 

-16.1* 

4 

2.5 

-3.1 

2.5 

-7.2 

5 

3.5 

0.8 

1.75 

-7.8 

6 

4.5 

4.1 

N/A 

N/A 

* This value wou 

d probably be changed by the integrator in time. 


Conclusion 

The analysis here shows that the basic PI controller has good stability margins. It typically 
produces a single stage change in a time interval of 0.25 to 0.35 seconds. It was shown 
that the dips in the fuel flow due to the speed of the stage change can be reduced by 
modifying the demanded fuel flow. The position errors for a stage change were generally 
less than ±0.010 inches. The fuel nozzle pressure drop stayed above the requirement of 20 
psi in all of the runs. However, simulations of a complete set of operating points is 
needed to better determine whether the pressure requirement will always be met and this 
is beyond the scope of this work. 

The only area of performance in which the basic controller did not perform well was with 
the uncommanded piston movements due to large changes in the manifold pressure 
during quiescent periods of the demanded fuel flow. For these uncommanded 
movements, a pressure feedback controller was proposed and tested with simulations. It 
was shown to be very effective in preventing these uncommanded movements. 


Issues for Further Study 

Before the proposed control is used, a dispersion analysis should be performed and a 
decision should be made to determine which method should be used to reduce the 
oscillations for stage six. The dispersion analysis needs to determine how much effect the 
spring rate and preload will have on the controller performance. The oscillations of stage 
six can be reduced by using a pressure transducer with a bandwidth higher than 8 hz or by 
lowering the gain for stage six. 
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It is also recommended to derive a transfer function model of the CPFN and CPFN PSV 
loop when the piston is static. This model could then be used to generate frequency 
domain plots and root locus plots for the static controller. These could in turn be used to 
quickly design and evaluate a pressure feedback controller that may outperform the one 
given here. 
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function [stab,gmdb,pm,wgain,wphase]=cpfnstab(kp,ki,plotopt,printopt) 

% [stab,gmdb,pm,wgain,wphase]=cpfnstab(kp,ki,plotopt,printopt) 

% 

% determines the stability margins of the CPFN loop 
% based on 2 A n parameter combinations 

% 

% parameters: PI, dp=(Ps-Pl), Qpsveq 

% 

% 

% inputs: 

% kp = proportional gain 

% ki = integral gain 

% plotopt= 0 no plot 

% 1 plot 

% printopt = 0 don't print plot file 

% 1 print plot file 

% 

% 

% outputs: 

% stab = stability flag, value = stable or unstable 

% gmdb = gain margin in decibels 

% pm = phase margin 

% wgain = frequency at which the gain margin is obtained (rad/sec) 

% wphase = frequency at which the phase margin is obtained (rad/sec) 
% 

% written by Rick Lalonde written 8-31-98 revised 9-10-98 

if (nargin < 1), kp=2200; end 
if (nargin < 2), ki=24; end 
if (nargin < 3), plotopt=0; end 
if (nargin < 4), printopt=0; end 

area=0.785; 

hl3=0.3935; hl4=0.18; T=0.01; 
fl9=0.6065; hl9=0.1804; il9=0.2131; 

w=logspace(- 1 ,2,300)'; 

% initialize "stab" to stable, if unstable roots are found this is changed 
stab-stable'; 


% 

plmin=40; plmax=200; 
dpmin=300; dpmax=500; 
eqmin=l; eqmax=2; 
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% 

case=0; 
for ii=l:2 
for jj=l:2 
for kk=l:2 
case=case+l; 

if (ii==l), pl=plmin; else pl=plmax; end 
if (jj=l), dp=dpmin; else dp=dpmax; end 
if (kk=l), eq=eqmin; else eq=eqmax; end 

ps=pl+dp; 

% t.f. from CPFN velocity to position 
nl=[0T];dl=[l-l]; 

% sensor t.f. 

n2=[0 hl9+il9]; d2=[l -fl9]; 
num=conv(nl,n2); den=conv(dl,d2); 


% 

if (eq=l), 

kpps v=sqrt( 1 00/(p 1 -20)) ; 
gainpsv=l/kppsv; 
end 

if (eq=2), 

kppsv=sqrt( 1 00/(ps-p 1 )); 
gainpsv=l/kppsv; 
end 

%kppsv=l; 

% controller t.f. 

n3=[kp*kppsv ki*T-kp*kppsv]; d3=[l -1]; 
num=conv(num,n3); den=conv(den,d3); 

% t.f. from TMC4 to Qpsv 
n4=gainpsv*[0 hl3*hl4]; d4=[l -(l-hl3)]; 
num=conv(num,n4); den=conv(den,d4); 

% gain from Qpsv to velocity 
n5=-l / (30*area); d5=l; 
num=n5*num; 

% calculate closed loop t.f. 
num=-num; 
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numcl=num; dencl=num+den; 


% calculate roots 
poles=roots(dencl); 

if any((abs(poles) >= 1)), stab='unstable'; end 

% plot Nichols chart 
[m,ph]=dbode(num,den,T,w); 
mdb=20*logl0(m); 

if (plotopt=l | printopt— 1), 
figure(case); 

plot(-180,0,'*',ph,mdb); grid; 

xlabel('PHASE (deg)'); ylabel('MAG (db)'); 

title([' "cpfnstab.m" pl= ’,num2str(pl),' ps= ’,num2str(ps), ... 

' eq= ',num2str(eq),’ case= ',num2str(case),' ',date]); 
end 

if (printopt=l), 

print -dps -append c:\cpfnstab.ps 
end 

[gm,pm(case, 1 ),wgain(case, l),wphase(case, 1 )]=margin(m,ph,w); 
gmdb(case,l)=20*logl0(gm); 

plvect(case,l)=pl; psvect(case,l)=ps; eqvect(case,l)=eq; 
kpvect(case, 1 )=kppsv; 

end % ii loop 
end % jj loop 
end % kk loop 

dispCpl ps eq KPpsv '); 

[plvect psvect eqvect kpvect] 

% 

% 
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Xcpfn demand and sensed col BV and Bx 




Appendix B - Position Control with Breakaways, Stages 1 to 2 
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Xcpfn demand and sensed, col BV and Bx 
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Figure C-l 

Appendix C - Position and Pressure Control, Stages 1 to 6 
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Controller switch & position error, col GO & BY 
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Figure C-2 

Appendix C - Position and Pressure Control, Stages 1 to 6 
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Controller switch & pressure error, col GO & GE 
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Figure C-3 


Appendix C - Position and Pressure Control, Stages 1 to 6 
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FMV demand, FMV flow & WF36, col B, O, & EP 



FMV demand - FMV add - WF36, col FJ 
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Figure C-4 

Appendix C - Position and Pressure Control, Stages 1 to 6 
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Xcpfn error col BY 
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Appendix C - Position and Pressure Control, Stages 1 to 6 
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Xcpfn demand and sensed, col BV and Bx 
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Figure D-l 

Appendix D - Position and Pressure Control With High Bandwidth Pressure Transducer, 

Stages 1 to 6 
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Figure D-2 


Appendix D - Position and Pressure Control With High Bandwidth Pressure Transducer, 

Stages 1 to 6 
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Figure D-3 


Appendix D - Position and Pressure Control With High Bandwidth Pressure Transducer, 

Stages 1 to 6 
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FMV demand, FMV flow & WF36, coi B, O, & EP 





Figure D-4 


Appendix D - Position and Pressure Control With High Bandwidth Pressure Transducer, 

Stages 1 to 6 
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Figure D-5 

Appendix D - Position and Pressure Control With High Bandwidth Pressure Transducer, 

Stages 1 to 6 
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Xcpfn demand and sensed, col BV and Bx 




Appendix E - Position Control With Decreasing Flow, Stages 6 to 1 
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Figure E-2 


Appendix E - Position Control With Decreasing Flow, Stages 6 to 1 
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Figure E-3 


Appendix E - Position Control With Decreasing Flow, Stages 6 to 1 
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FMV demand, FMV flow & WF36, col B, O, & EP 




Figure E-4 

Appendix E - Position Control With Decreasing Flow, Stages 6 to 1 
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Xcpfn error col BY 
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Figure E-5 


Appendix E - Position Control With Decreasing Flow, Stages 6 to 1 
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Xcpfn demand and sensed, col BV and Bx 
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Figure F-l 

Appendix F - Position Control With fuel Demand Modification, Stages 2 to 4 
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Figure F-2 


Appendix F - Position Control With fuel Demand Modification, Stages 2 to 4 
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Figure F-3 

Appendix F - Position Control With fuel Demand Modification, Stages 2 to 4 
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FMV demand, FMV flow & WF36, col B, O, & EP 
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Figure F-4 

Appendix F - Position Control With fuel Demand Modification, Stages 2 to 4 
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Figure F-5 

Appendix F - Position Control With fuel Demand Modification, Stages 2 to 4 
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Attachment III: Neural Networks for NOx Emission 
Prediction - Final Report 
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Introduction 


The purpose of this work was to evaluate the efficacy of artificial neural networks (NN) 
for prediction of oxides of nitrogen produced by aeronautical or aeronautically derived 
turbine engines. The terminology “oxides of nitrogen” covers both nitric oxide (NO) and 
nitrous oxide (NO 2 ) and is denoted by the formula NO x where x is a variable. The 
usefulness of the NNs should be determined mostly by their prediction error, however, 
complexity and ease of prediction are also important. The prediction of the NNs is 
compared here to the prediction errors of the previous prediction standard: empirically 
derived exponential models. 

The data available for this phase of the work is from two engine tests of two different 
LM6000 engines fueled by natural gas. In each data set there are 36 inputs and 6 outputs. 
These are listed in Appendix A. The inputs are mostly temperatures, pressures, and flows. 
The outputs are the raw and adjusted values of NO x , carbon monoxide, and unbumed 
hydrocarbons. The only output considered for this work is NOX15, the adjusted NO x 
value. 

Goals of this work were to 

1 Evaluate different NN structures and training methods and determine a structure 
and training method that will work well for predicting NOX emissions. 

2 Determine how many neurons are needed to get good prediction results. 

3 Determine how large the data set needs to be to train the network effectively. 

4 Make the prediction error small. 

The main goal is to generate one model that is good for all LM6000 engines that meet the 
design specifications. Therefore, the model should perform well on LM6000 engines that 
it has not been trained with. Good performance only on the engine used for training is not 
considered sufficient. 

NO x Measurement Error and Surface Roughness 

The NN can not be expected to produce a model with an error smaller than the error of 
the output that was given to it during training. For this reason, it was desired to estimate 
the sensor measurement error. The error was estimated in the following manner. 

The inputs were normalized and for each operating point in the data set, the point with 
the closest input vector was determined. A plot was then made of the change in the value 
of the output as a function of the distance between the input vectors. If two inputs exist 
that produce different outputs, then a pathology exists and maybe caused by one of the 
following: the input set does not contain all of the necessary inputs (that affect the 
output), or the actual outputs are the same and the sensor error causes the sensed outputs 
to be different. For data set 1 12, there were no pathological points when all of the inputs 
were used, see Figure B-l. In fact, in Figure B-2, a subset of six inputs is shown not to 
have any pathological points. For lack of a better name, these plots will be referred to 
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here as roughness plots as the variation of the points at a particular distance show the 
roughness of the output surface over the input. 

When two straight lines are drawn to represent the top and bottom of the envelope of 
points on the roughness diagram, they become closer together as they move from right to 
left. The distance between the lines at the left is an estimate of the sensor measurement 
error. For data set 1 12 this error seems to be conservatively less than 0.5 PPM for NO x . 

Theoretical Equations for NO x Production 

According to [1], oxides of nitrogen can be produced in gas turbine engines by at least 
four different mechanisms: Thermal NO, the Nitrous Oxide Mechanism, Prompt NO, and 
Fuel NO. Balanced chemical equations are listed in that reference. However, the amount 
of production of the oxides of nitrogen by these four mechanisms depends on the fuel, the 
type of fuel atomization, the droplet size, temperatures at various stages, pressures at 
various stages, the equivalence ratio, the residence time, etc. Because of all of these 
factors, the determination of a set of equations for the amount of NO x produced is 
currently, nearly intractable, even for a specific engine. To date, researchers have 
modeled the amount of NO x with equations using a subset of the above factors. The 
constants in the equations have been determined by curve fitting and regression 
techniques. This paper considers two classes of models. The first class consists of models 
in which the NO x is an exponential function of temperatures. The second class of models 
is artificial neural network models. 

Single Input Exponential NO x Model 

NO x is plotted versus TFLPDF for data set 1 12 in Figure D-l. (This plot contains only 
the points for BRNDMD=15.) Since the data for a single combustion chamber is 
expected to be exponential, this plot was expected to have an exponential relationship 
with additive “noise” due to, at the least, the inner and outer flame temperatures. By 
inspection an exponential curve can be drawn through some of the data starting at about 
(3100,16) and ending at (3450,60). Points were chosen along this path every 50 degrees 
and then were curve fit to the equation 

NO x = e (a * temp + b) 


The fit resulted in the values of a=0.004016 and b=-9.7333. 

Three Input Exponential NO x Model 

General Electric’s Aircraft Engine (GEAE) Division has a NO x prediction model for the 
LM6000 engine. The model is 


NO x s (1 / sum(wp(i))) Z wp(i) e (cp(,) * (tfl(i) + dtflp(l » +mup(i)) 
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where: 

tfl(l) = outer flame temperature 
tfl(2) = pilot flame temperature 
tfl(3) = inner flame temperature 

wp(i) = flow fraction for outer, pilot, and inner, respectively 

Regularization and Enhancing Performance 

Neural nets are generally thought to be very capable of producing targets close to those 
of the training set. (The set of untrained inputs and outputs is often referred to as the 
complimentary set.) Unfortunately, this is no guarantee that the nets will produce good 
outputs in response to inputs they have not seen. This can be attributed to several causes. 
One is that the training set is not sufficiently large, so that the inputs in the 
complimentary set are more like extrapolations of the training set than interpolations. The 
second reason also has to do with training sets that are not large enough. If the number of 
parameters of the net is roughly the same size or larger than the number of inputs, then 
the network can “overfit”. Sometimes overfitting is so apparent that the error on the 
training set is practically zero, while the error on the complimentary set is very large 
because the model has severe fluctuations between the points in the training set. 

A net is said to regularize well if the errors on the complimentary set are not much larger 
than those of the training set. Therefore, regularization is a property that is desired of 
most networks. 

For this particular application (gas turbine engines), the training set could be data from 
one engine and the complimentary set could be data from the same engine. In this case, 
training of the NN could happen before the unit is sold. The training set could also be 
data from a test or reference engine and the complimentary set could be data from other 
engines. For either of these cases, the most important set to perform well on is the 
complimentary set. That is why this analysis focused heavily on reducing the prediction 
error of the complimentary set by reducing overfitting and by providing good 
generalization. Making the prediction error of the training set small, by itself, is of little 
value. 

Early stopping and using a validation set are two ways to avoid or reduce overfitting. 
They will be described briefly below, followed by two techniques and a network 
architecture that were developed to provide good generalization. There are many papers 
in the literature on techniques for regularization and to avoid overfitting. The ones 
presented here are a subset of those that the author used for this task. 

Early Stopping 

The most simplistic approach to avoid overfitting, other than to select a small number of 
neurons, is to train the network for a small number of epochs, (training for one epoch 
means to present the training set to the network only one time). The hope is that with just 
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a few epochs, the net will learn the input/output relationship relatively well, but will not 
fit the noise in the training set. 

Validation 


Another common practice to avoid overfitting is to split the data into three sets: a training 
set, a validation, and a complimentary set. As the network is being trained, calculate the 
prediction errors for both the training and complimentary sets. The error of the training 
set will often be monotonically decreasing. When the error of the validation set starts to 
increase, then save the network parameters from when the validation set had the 
minimum error. The network is then tested on the complimentary set as usual. 

Training with Bayesian Regularization 

Another way to avoid overtraining and improve generalization is Bayesian Regularization 
(BR). This training algorithm is a modified version of the Levenberg-Marquardt 
algorithm in which the error function has an added term which is the sum of the squares 
of the network parameters. This added term places a penalty on the parameters as well as 
the prediction error. Thus as the network is being trained there is a trade-off between 
error of the training set and the size of the parameters. This reduces the effective number 
of parameters (i.e. effective size) of the network, thereby reducing the probability of 
overfitting. 

Radial Basis Networks 


There are two main types of radial basis networks: exact and inexact. The exact networks 
have exactly zero error on the training set because they have as many neurons in the input 
layer as data points in the training set. This is a severe drawback, but they generalize well 
do the architecture of their first layer. In the first layer, the input vector is compared to a 
vector stored in each neuron (these vectors comprise the radial basis). The distance 
between the input vector and the stored vector is computed. If the distance is zero, the 
output of the neuron is one, otherwise it is positive and less than one. Therefore, the 
neurons in the first layer only respond well to input vectors close to their stored vector. 
The second layer takes a weighted average of the targets of the first layer by using the 
outputs of the first layer as weights. 

The inexact radial basis networks have the same architecture as the exact radial basis 
nets, but fewer neurons. This causes them to have non-zero error on the training set, but 
they generalize well again because of their first layer. 

Input Jitter 

Another technique to produce a net that will generalize well is to add jitter to the inputs. 
This is done by taking each input of the training set and making “n” copies of the input 
each with additive (and independent) noise. For each artificial input use the exact same 
output that the real input had. The new training set will be (n+1) times the size of the 
original set. This will force the network to be relatively smooth around each of the inputs 
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in the training set, thereby providing better regularization. If several of the artificial 
inputs are very close together, then the network will generally produce an output for that 
region of input space that is close to the average of the outputs in that region. Since it is a 
goal to have many more data points than parameters in the net, these artificial points are 
often sufficient to meet the goal. 

Generalization Techniques Applied to the LM6000 Engine Data 

All five techniques described in the previous section to reduce overfitting and improve 
generalization were used in moderately thorough tests with the available engine data. 

The best results were obtained using input jitter and the Bayesian Regularization training 
method together. Therefore, only results obtained by using these two methods will be 
presented in this report. 

Calculated Pilot Ring Flame Temperature for Data Set 112 

The data set for engine serial number 1 12 has four burner demands: 0, 15, 25, and 40. 
Figure D-l shows the adjusted NO x , NOX15, for a burner demand of 15. The raw data is 
plotted with asterisks. It is seen that much of the data can be fit reasonably well to a 
curved line. The curved line in the plot is NO x = e (a * temp + b) with a=0.004016 and 
b= -0.73333. Since the points in the upper left portion of the graph are not approximated 
well by the curved line, the validity of the calculation for the pilot temperature became a 
small concern. The calculated pilot flame temperature is not a measurement, but is an 
estimation, like the other calculated flame temperatures. They are based on the demanded 
inner and outer ring flame temperatures and other variables. Therefore the relationship 
between NOX15 and the pilot temperature is only as good as the calculation of the pilot 
temperature. 

Figure D-2 shows the demanded and calculated inner ring flame temperatures for burner 
demands of 15 and 40. For burner demands of 15, it is seen that some of the predicted 
data points are far from the demand. This does not necessarily mean that the prediction is 
inaccurate. (Since the actual flame temperatures can not be measured in a practical 
manner, the accuracy of the calculation can only be inferred here and can not be 
determined.) To see the effect of this difference between the demanded and calculated 
data points, the calculated temperatures above the line in Figure D-2 were adjusted by 
120 deg F, creating plots D-3 and D-4. Figure D-3 shows the demanded and calculated 
inner ring flame temperatures after this modification. Figure D-4 shows the values of 
NOX15 versus calculated pilot ring flame temperatures with the modification to them. 
They were shifted 120 deg in the opposite direction, thereby keeping the average (and 
sum) of the flame temperatures the same before and after modification. The plot of 
NOX15 versus pilot flame temperature is much more like an exponential curve now. 

The relationship between NOX15 and the pilot temperature seems to be exponential 
with corrections for factors other than the pilot temperature. 
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Results for Data Set 112 


As a baseline for determining the useful of the NNs, the standard deviations (SDs) and 
means of the errors for the exponential models were calculated for Data Set 1 12. These 
are listed in Tables 1 and 2 for the raw data set and in Tables 3 and 4 for the data set with 
the previously mentioned modification of the calculated flame temperatures. The 
parameters for the three input exponential model were based only on a burner demand of 
40. This is why the table does not have entries for other burner demands for this model. 
The column labeled “ALL” designates all four values of the burner demand. It is seen 
that the modification of the calculated flame temperatures, for the burner demand value 
of 15, reduces the SDs of the error of both exponential models. This lends credibility to 
the modification. 

Results using a feed-forward NN with input jitter and the Bayesian Regularization 
training method are shown in Appendix E. Figure E-l shows the SDs of the error of the 
training and complimentary sets versus training size. To make this graph. Modified Data 
Set 112, which has 179 operating points, was randomly ordered. Then the first 20 points 
were used as a training set and the last 159 as the complimentary set. After this, the first 
30 were used as the training set. Therefore, 20 of the operating points were the same as 
the previous run. This was then completed until 160 points were used for the training set, 
thereby leaving only 19 points in the complimentary set. It is thought that fewer than 19 
points is not statistically significant. Figures E-2 through E-5 show the results with 
different initial random orderings. Finally, Figure E-6 is average of the previous five 
graphs. The mean and standard deviation for the largest training set size, namely 160, is 
given at the top of the plot and in Table 2. The SD of the training set seems to be settling 
around 2.4, but the SD of the complimentary set still seems to be decreasing and appears 
like it would go further below 3.0 if the data set were large enough The standard 
deviation of the neural network is 24% of the three input model and 55% of the one input 
model. 

Forty epochs were used for all of the results. Since jitter increased the number of data 
points by a factor of five, the ratio of data points to NN parameters was five times higher 
and over training was not a problem. 


Table 1 Standard deviations of the error for Data Set 112 without the modification. 


BURNER DEMAND 

0 

15 

25 

40 

ALL 

Single Input Exp. Model 

3.66 

8.14 

4.53 

5.36 

6.79 

Three Input Exp. Model 




4.98 

11.65* 

NN Model using Jitter 
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Table 2 Standard deviations of the error for Data Set 112 with the modification. 


BURNER DEMAND 

0 

15 

25 

40 

ALL 

Single Input Exp. Model 

3.66 

3.14 

4.53 

5.36 

4.67 | 

Three Input Exp. Model 




4.98 

10.67* 

NN Model using Jitter 





2 - 59 1 


Table 3 Means of the error for Data Set 1 12 without the modification. 


BURNER DEMAND 

0 

15 

25 

40 

ALL 

Single Input Exp. Model 

-2.59 

-4.47 

3.04 

0.10 

-1.27 

Three Input Exp. Model 




2.42 

4.67* 

NN Model using Jitter 







Table 4 Means of the error for Data Set 1 12 with the modification. 


BURNER DEMAND 

0 

15 

25 

40 

ALL I 

Single Input Exp. Model 

-2.59 

0.49 

3.04 

0.10 

0.42 

Three Input Exp. Model 




2.42 

5.18* 

NN Model using Jitter 





-0.35 


Results for Data Set 118 

Figures F-l to F-5 show the results of five different initial random orderings of Data Set 
118. Figure F-6 gives the average of the previous five plots. In these figures, it is seen 
that the SDs of the training set are smaller than those of Data Set 112, even for the same 
training set size. This is also true for the exponential models as can be seen by Tables 5 
and 6. The reason for this is not known. The SD of the NN models is 57% of the three 
input exponential model 29% of the one input exponential model. These percentages are 
quite similar to those of Data Set 1 12, which are 55% and 24% respectively. 


Table 5 Standard deviations of the error for Data Set 1 1 8. 


BURNER DEMAND 

0 

15 

25 

40 

ALL I 

Single Input Exp. Model 

3.12 

1.86 

5.50 

2.35 

3.62 

Three Input Exp. Model 




1.13 

7.19* 

NN Model using Jitter 





2.05 


Table 6 Means of the error for Data Set 1 1 8. 


BURNER DEMAND 

0 

15 

25 

40 

ALL 

Single Input Exp. Model 

-1.36 

-3.32 

1.33 

-2.47 

-1.45 

Three Input Exp. Model 




3.43 

-0.26* 

NN Model using Jitter 





0.06 


* The three input exponential model constants were only intended for a burner demand of 40. 
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Results for Data Sets 112 & 118 Combined 


Data Sets 112 and 118 were combined in two different ways for this work. The first way 
was to randomly order each set and then combine the sets to make the first 179 operating 
points of the new set to be those of 1 12 and to make the last 66 points to be those of 1 1 8. 
The second way is to merge the two sets and then to randomly order the combined set. 
The first way should present a more difficult test to the NN, because for small training set 
sizes, the training set will not have any of Data Set 1 18 in it, but the complimentary set 
will. In other words, the second way is easier because, the NN is trained on both sets and 
tested on both sets. 

The results from combining the sets the first way is shown in Figures G-l to G-6. The 
resulting mean is 0.03 and the SD is 2.77. Figure G-6 is quite important. It shows that 
when the training set size is 170 and the NN has been trained only on Data Set 112, that it 
performs as well on Data Set 118 as a NN trained only on Data Set 118 (compare to 
complementary set of Figure F-6). This shows that the underlying relationship between 
the inputs and outputs is quite similar with both LM6000 data sets, as it should be. The 
network is seen to generalize well; the SDs of the complimentary set are not much higher 
than those of the training set. The input set for all of the results presented in this paper is 
the set of variables [9 15 16 17 20 22] from page A-l. 

The results from combining the data sets the second way are shown in Appendix H. The 
resulting mean is 0.19 and the SD is 2.65. Even though theoretically, the first way of 
combining the sets should pose a more difficult test, the results are relatively close. One 
has a lower mean, the other has lower SD. The SDs given here from the combined set are 
similar to the SDs for Data Set 112. 

Table 7 compares the single input exponential model to the NN models. The three input 
exponential model of [2] is not compared here because it was only intended for one value 
of the burner demand and it has different coefficients for each engine. Therefore, it is 
much more limited. 
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Table 7 Results from combining Data Sets 1 12 & 

L18. 

MODEL 

MEAN 

SD 

SD 

RATIO 

Single Input Exp. Model 

-0.08 

4.39 

100% 

NN randomization & merging of the data sets (App. G) 

0.03 

2.77 

63% 

NN merging & randomization of the data sets (App. H) 

0.19 

2.65 

60% 


It is seen from the table, that the SDs of the errors of the NN model are roughly 60% of 
those of the single input exponential model. 

Conclusion 

The most important contribution of this work was to determine how to generate a NN that 
would generalize well, so that it could be trained on one engine and perform well on 
another. To keep this report concise, results from the different methods (listed in this 
report) to provide generalization and to enhance performance were not given here. 
Relatively thorough tests were conducted using the previously described methods: early 
stopping, validation, radial basis networks, training with Bayesian Regularization, and 
input jitter. It was determined that the last two work the best for this application and 
therefore only those results were shown. 

This work has shown that by using input jitter and by training feed-forward networks of 
only 12 neurons using Bayesian Regularization, that NO x prediction errors of about 3.0 
can be achieved for data sets that the NN has not been trained on. Furthermore, due to the 
use of input jitter and the BR training method, the relationships between the inputs and 
the output are smooth. This greatly reduces the probability of the NN providing a very 
poor NO x prediction on any input that can be considered as an interpolation of the 
training inputs. Since jitter increased the number of data points by a factor of five, the 
ratio of data points to NN parameters was five times higher and overfitting was not a 
problem. 

It is always preferably to have many data sets to draw inferences from. However, it 
appears that if the two used for this study are typical, that a NN trained on one LM6000 
engine (with 180 operating points or more) can provide good performance on other 
LM6000 engines when input jitter is used and the network is trained with Bayesian 
Regularization. 
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Appendix A - Variable Description for Data Sets 112 & 118 


1 

esn 

engine serial number 


2 

rdg 

Data point reading number 


3 

tdate 

Date 


4 

bmreq 

Burner configuration SS target 

n/d 

5 

bmdmd 

Burner configuration demand 

n/d 

6 

nggsel 

Gas Generator speed 

rpm 

7 

nggdot 

Gas generator speed rate of change 

rpm/sec 

8 

t2sel 

Gas generator inlet temperature 

degF 

9 

t3est 

Compressor Discharge Temperature 

degF 

10 

t54sel 

Turbine dischage temperature 

degF 

11 

tflcycs 

Smoothed calculated bulk flame temperature 

degF 

12 

tfldmd 

Bulk flame temperature demand 

degF 

13 

tflirf 

Inner ring flame temperature demand 

degF 

14 

tflorf 

Outer ring flame temperature demand 

degF 

15 

tflidf 

Calculated inner ringflame temperature 

degF 

16 

tflpdf 

Calculated pilot ring flame temperature 

degF 

17 

tflodf 

Calculated outer ring flame temperature 

degF 

18 

tflpct 

Bulk flame percent (min=0, max=100) 

% 

19 

p2sel 

Gas generator inlet pressure 

psia 

20 

ps3sel 

Compressor dischage pressure 

psia 

21 

p54sel 

Turbine Discharge pressure 

psia 

22 

vvf36dmd 

Fuel flow demand 

PPh 

23 

wflabc 

Inner fuel flow reference per staging valve 

PPh 

24 

wfoabc 

Outer fuel floe reference per staging valve 

PPh 

25 

wfinr 

Calculated inner ring fuel flow 

PPh 

26 

wfpil 

Calculated pilot ring fuel flow 

PPh 

27 

wfotr 

Calculated outer ring fuel flow 

PPh 

28 

wfinrs 

Smoothed calculated inner ring fuel flow 

PPh 

29 

wfpils 

Smoothed calculated pilot ring fuel flow 

PPh 

30 

wfotrs 

Smoothed calculated outer ring fuel flow 

PPh 

31 

fwa36 

Tflame airflow at station 36 

pps 

32 

wacfli 

Inner airflow burned in combustion 

pps 

33 

wacflo 

Outer airflow burned in combustion 

pps 

34 

wacflp 

Pilot airflow burned in combustion 

pps 

35 

wb26q 

Eight stage fractional bleed flow demand 

n/d 

36 

wb3q 

CDP fractional bleed flow demand 

n/d 

37 

dwb36 

Total bleed flow demand reflected to sta 36 

pps 

38 

st8sel 

Eight stage bleed valve sensed position 

% stroke 

39 

cdpsel 

CDP bleed valve sensed position 

% stroke 

40 

nox 

Raw Nox level 

ppm 

41 

noxl5 

Nox level corrected to 15% 02 

ppm 

42 

CO 

Raw Co level 

ppm 

43 

col5 

Co level corrected to 15% 02 

ppm 

44 

uhc 

Raw unburned hydrocarbonlevel 

ppm 

45 

uhcl5 

Unbumed hyhdorcarbon level corrected to 15% 02 

ppm 

46 

lhvsel 

Fuel lower heating value 

btu/lbm 

47 

sgsel 

Fuel specific gravity 

n/d 
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Roughness Plot for data set 1 12 with all inputs 



Figure B-l Roughness plot with all of the inputs. 


Roughness Plot for data set 112 with inputs [9 15 16 17 20 22] 



distance between normalized input vectors 


Figure B-2 Roughness plot with inputs 9, 15, 16, 17, 20, & 22 . 
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Appendix C - Parameters for Three Input Exponential NOX Model 


General Electric determined parameters for the three input exponential model for data 
sets 1 12 and 118. The Parameters are intended only for a burner demand (BRNDMD) of 
40. The model equation is repeated below: 

NO x = (1 / sum(wp(i))) I wp(i) e (c P (,) * (tfl(1) + dtfl P (,)) +mu P (,)) 

The parameters are: 


dtflp 112 = [1.0795e2 

-2.201 lei 

6.5252el] 

dtflp 118 = [9.5405el 

-3.5124el 

-9.443 lei] 

cpll2 = [1.2609e-2 

5.6860e-3 

5.6840e-3] 

cpll8 = [4.5500e-3 

5.6860e-3 

5.6840e-3] 

mupll2 = [-4.2960el 

-1.6926el 

-1.6642el] 

mupll8 = [-1.2581el 

-1.6926el 

-1.6642el] 


wpll2 = wpll8 = [0.295 0.320 0.139] 
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Figure D-2 Inner ring flame temperatures for Data Set 1 12 without the modification: 
calculated versus demanded temperature 
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Burner Demand = 15 




Figure D-3 Inner ring flame temperatures for Data Set 1 12 with the modification: 
calculated versus demanded temperature 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 



Training Set Size 


Figure E-l Results of random ordering #1 of Data Set 1 12 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 



20 40 60 80 100 120 140 160 

Training Set Size 


Figure E-2 Results of random ordering #2 of Data Set 1 12 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 
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Figure E-3 Results of random ordering #3 of Data Set 1 12 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 



Training Set Size 


Figure E-4 Results of random ordering #4 of Data Set 1 12 
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Std. Dev. of the Error 


loopjrt.m * = train o = comp. 



Training Set Size 


Figure E-5 Results of random ordering #5 of Data Set 1 12 
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Std. Dev. of the Error 


Data Set 112 random average of 5 runs with jitter std. dev.= 2.5855 mean= -0.34984 



Figure E-6 Results of 5 random orderings of Data Set 1 12 
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Std. Dev. of the Error 



Training Set Size 


Figure F-l Results of random ordering #1 of Data Set 118 
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Dev. 


loopjit.m * = train o = comp. 



Training Set Size 


Figure F-2 Results of random ordering #2 of Data Set 1 18 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 



Training Set Size 


Figure F-3 Results of random ordering #3 of Data Set 1 1 8 
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Std. Dev. of the Error 


loopjit.m * = train o = comp. 
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Figure F-4 Results of random ordering #4 of Data Set 1 18 
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Std. Dev. of the Error 


5 


Data Set 118 random average of 5 runs with jitter std. dev. = 2.0457 mean= 0.056381 



Figure F-6 Results of 5 random orderings of Data Set 1 18 
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Training Set Size 


Figure G-3 Results of random ordering #3 of the separate Data Sets 1 12 & 1 1 8 
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Training Set Size 


Figure G-4 Results of random ordering #4 of the separate Data Sets 112 & 1 18 
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Std. Dev. of the Error 



Training Set Size 


Figure G-5 Results of random ordering #5 of the separate Data Sets 1 12 & 118 
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Std. Dev. of the Error 


Data Sets 112 & 118 separate, average of 5 runs with jitter std. dev.= 2.7731 mean= 0.025043 



Figure G-6 Results of 5 random orderings of the separate Data Sets 1 12 & 118 
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Std. Dev. of the Error 
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Figure H-2 Results of random ordering #2 of the combined Data Set 1 12 & 1 18 
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Std. Dev. of the Error 
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Figure H-3 Results of random ordering #3 of the combined Data Set 1 12 & 1 18 
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Std. Dev. of the Error 



Training Set Size 


Figure H-6 Results of 5 random orderings of the combined Data Set 1 12 & 1 18 
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Std. Dev. of the Error 



Training Set Size 


Figure H-7 Results of 1 1 random orderings of the combined Data Set 1 12 & 1 18 
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Table J-l Matlab functions written for this neural network task. 


DATE 

*.MFILE 

DESCRIPTION 

10-2-98 

nnscater 

output scatter plots with different colors and symbols 

10-2-98 

datacorr 

correlation of vectors 

12-1-98 


exponential model 

12-1-98 

identexp 

identifies an exponential model 

12-2-98 

datadist 

calculates distance apart of input vectors 

12-3-98 

datamodl 

interpolate NN input data (input data conditioning) 

12-3-98 

sdcomp 

calculates mean and SD of complimentary set 

12-4-98 

nettest 

short (one page) & simple NN setup and training function 

12-10-98 


divides three temps by flows 

12-18-98 

inscater 

plots inputs with a color scatter 

12-18-98 


implements “early stopping” 

12-21-98 

trainsiz 

plots SD of error versus training set size 

1-6-99 

netsize 

plots SD of error versus NN size 

1-22-99 


averages “n” NN models 

1-26-99 

radbtest 

short and simple implementation of radial basis networks 

1-26-99 

sizegmn 

short and simple implementation of gen. regression NNs 1 

1-27-99 


makes “roughness” plots 1 

1-28-99 

findexp 

finds operating points that are near single input exp. model 1 

2-1-99 

regress 

regression jj 

2-2-99 


GE’s three input exponential model 1 

2-3-99 

bee— i 

adds jitter to the input data set 

2-8-99 

H9HI 

calls jitter and nettest repeatedly to make plot versus size 

2-11-99 

sortdata 

sorts based on burner demand and plots 

2-11-99 

pilotemp 

modifies Data Set 1 12 pilot temperature f 


The most important functions are in boldface type. 
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Attachment IV: FMUX Dynamic Modeling - Summary 

Memorandum 
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SUBJECT : Addition of FMU to Active Combustor Control Model 
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PROJECT NUMBER: X83753 
FROM: Daniel McLevige 
DATE: August 28, 1998 
TO: Bill Myers, GEAE 
cc: Ted Hommema 

1.0 ABSTRACT 

The FMU and pump were added to the Active Combustor Control Simulink dynamic model to 
check for interactions between the two systems. No significant dynamic interactions were found, 
but a steady state total flow change of up to 1.5% was found due to bypass valve droop at 125% 
flow split. 

2.0 CONCLUSIONS & RECOMMENDATIONS 

a) The FMU and pump were added to the Active Combustor Control Simulink dynamic model and 
runs made stepping both total flow and flow split. 

b) Flow split step behavior was essentially identical to previous runs of the model before the FMU 
was added. 

c) No adverse dynamic interaction between the HMUX flow splitting function and FMU flow 
metering function was seen. 

d) As flow was uptrimmed on one nozzle, the resultant increase in back pressure caused total flow 
drop because of bypass valve droop. Up to 1.5% total flow drop was seen at 125% nozzle flow 
uptrim. If this is deemed a concern, a closer look using a steady state model should be done, since 
the dynamic model is not optimized for precisely predicting steady state performance. 

e) Since the simulation does not model all 30 nozzles individually, it was found to give improper 
response of flow split to total flow steps. To model all 30 nozzles in detail is beyond the capability 
of the hardware/software used. 

3.0 DISCUSSION 

The Simulink model of the active combustor control HMUX has been used to predict dynamic 
behavior. Model runs involving step response in requested flow split between nozzles have been 
presented at the October 9, 1997 and May 27, 1998 design reviews at GEAE - Evendale. 
Representative model run plots were included in the hand-out books for the two reviews. At the 
last review, GEAE requested that the FMU be added to the model to look for interaction between 
the flow splitting function of the HMUX and the fuel metering function of the FMU. This was 
done and selected runs of the combined model made. 


NASA/CR— 2006-2 14244 


341 



Woodward Governor Company 
Rockford, Illinois 

3. 1 Model Summary 


EN831562 


Components included in the FMU model are as follows (see schematic Figure 1): 

a) Bypass valve - mass, damping, spring rate, damping orifices, nonlinear flow equation 

and flow force equation at all ports 

b) Pressurizing valve - mass, damping, spring rate, friction, damping orifice, non linear 

flow equation and flow force equation 

c) Boost pump - pressure as function of speed and flow 

d) Gear pump - flow as function of speed with internal leakage proportional to pressure 

rise 

e) FMU internal leakage - proportional to square root of Ps-Pb 

f) Fuel valve - port flow equation only 

g) PI, P2, and P22 volumes compressibility 

The active combustor control HMUX model was essentially the same as described in the previous 
design reviews, and included the following components (see schematic Figure 2): 

a) One HMUX only 

b) Force motor & first stage pilot valve 

c) Rotating sleeve (2 ports in detail, rest simplified) 

d) Hydraulic turbine - speed follows Ps-Pb with a time constant 

e) Pressure transducer 

f) Throttle valves (2 in detail, other 28 lumped together and simplified) 

g) Nozzles (2 in detail, other 28 lumped together and simplified, with option for turning 

nozzles off at low flows to maintain minimum nozzle delta P) 

h) Pc regulator - not used for runs made for this report (Assumed Ps supply to impulse 

turbine and HMUX first and seconds stages) 

i) Digital controller (2 nozzles plus delta P in detail, other 28 nozzles lumped together) 

j) Fluid compressibility and inertia in selected locations 

3.2 Model Runs 

Model runs were done at three total flows: 4,500 pph, 21,000 pph, and 29,000 pph. Similar cases 
were run as had previously been done with the model without FMU. One minor difference was 
that for the combined model runs, nozzles were shut off to maintain 30 psid nozzle delta P for the 
75% flow split case. The previous work had set the minimum nozzle delta P at 30 psid for the 
nominal 100% flow split case, meaning that delta P became excessively low for the downtrim 
cases. This change affected the 4,500 pph point only. 

Three cases were looked at for each flow level: 10% p-p flow split steps, + 1 - 25 % large flow split 
steps, and 10% p-p total flow steps. The flow split steps were compared to earlier model runs 
without the FMU and found to be essentially identical for the 21,000 pph and 29,000 pph cases. 
The 4,500 pph case was a little different because of the difference in the nozzle delta P. In all 
cases, no adverse dynamic interaction between the FMU flow metering and HMUX flow splitting 
functions was seen. 

However, a steady state effect was seen. As flow split changed, it changed the back pressure on 
the bypass valve for the cases off the pressurizing valve. Since the bypass valve does not contain 
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an integrator, this causes some delta P droop and flow error. While the model is optimized for 
showing dynamic rather than steady state performance, the bypass valve model is detailed enough 
that it should be close for steady state also. 

Model runs both with and without HMUX Pc regulator were done, but for brevity only plots from 
runs without Pc regulator are presented in this report. The runs with Pc regulator were similar. 

Plots of representative runs are attached as follows: 

a) Figure 3: 10% Flow Split Steps at 21,000 pph; Nozzles 1 & 2 Opposite Directions 

This shows a split time constant of about 3/4 seconds causing about 0.3% p-p total flow 
disturbances. 

b) Figure 4: +/- 25% Full Range Flow Split Steps at 21,000 pph. Nozzles 1 & 2 Opposite 
Directions 

This is a very unrealistic case, and is run only to prove that the system is stable with any 
conceivable step size. It shows total flow effects of 1.5% p-p. This is a steady state droop 
effect caused by the increased back pressure affecting the proportional bypass valve. 

c) Figure 5: 10% p-p Total Flow Steps from 21,000 pph 

This shows the total flow steps causing about a 1.5% temporary error in flow split. This 
is believed to be a model anomaly that is not expected to occur in the real system. In the 
model, only two nozzle channels were modeled in detail; the other 28 were lumped together 
and simplified. In the real system, all nozzle channels are identical, and a common mode 
effect such as stepping total flow should affect all equally. Therefore, for the case of all 
channels at 100% split, no significant change in flow split should occur during a total flow 
transient. 

d) Figure 6: 10% p-p Total Flow Steps from 4,500 pph 

This shows the total flow steps causing about a 3% temporary error in flow split. Again, 
this is believed to be an artifact of the model and should not happen in the real system. 

e) Figure 7: 10% p-p Total Flow' Steps from 29,000 pph 

This shows the same effect, in this case about a 4% temporary flow split error, which 
again should not be real. 
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Figure 1 HMUX Model Schematic 
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Figure 3 Flow Split Steps at 21,000 pph 
1 0% Steps Nozzles 1 & 2 Opposite Directions 
Variables vs. time: 

1) Percent flow split - channel # 1 

2) Percent flow split - channel #2 

3) Percent flow split - channels #3-30 lumped together 

4) Total flow at pressurizing valve discharge - percent 
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Figure 4 Full Range Flow Split Steps at 21.000 pph 
+/- 25% Steps, Nozzles 1 & 2 Opposite Directions 
Variables vs. time: 

1) Percent flow split - nozzle #1 

2) Percent flow split - nozzle #2 

3) Percent flow split - nozzles #3-30 lumped together 

4) Total flow at pressurizing valve discharge - percent 
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Figure 5 Total Flow Steps from 21,000 pph 
10% p-p Steps 
Variables vs. time: 

1 ) Total flow at pressurizing valve discharge - percent 

2) Percent flow split - nozzle # 1 

3) Percent flow split - nozzle #2 

4) Percent flow split - nozzles #3-30 lumped together 
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Figure 6 Total Flow Steps from 4,500 pph 
10% p-p Steps 
Variables vs. time: 


1) Total flow at pressurizing valve discharge - percent 

2) Percent flow split - nozzle #1 

3) Percent flow split - nozzle #2 

4) Percent flow split - nozzles #3-14 lumped together (other 16 nozzles closed off) 
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Figure 7 Total Flow Steps from 29,000 pph 
10% p-p Steps 
Variables vs. time: 

1) Total flow at pressurizing valve discharge - percent 

2) Percent flow split - nozzle # 1 

3) Percent flow split - nozzle #2 

4) Percent flow split - nozzles #3-30 lumped together 
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